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FUNDAMENTAL  STUDIES  ON  EROSION  IN  MPD  THRUSTERS 


Grant  No.  AFOSR-87-0360 

Final  Technical  Report 

V.  V.  Subramaniam 

Department  of  Mechanical  Engineering 
The  Ohio  State  University 
Columbus,  Ohio  43210 


1.0  INTRODUCTION 

The  purpose  of  this  research  is  to  understand  and  quantify  the  mechanisms 
responsible  for  evaporative  erosion  in  steady  state  magnetoplasmadynamic  (MPD) 
thrusters.  This  is  necessary  in  order  to  predict  thruster  characteristics  and  lifetimes  for 
a  given  design.  The  Back-EMF  theory  of  Onset  has  been  refined  and  expanded  to 
enable  prediction  of  erosion  rates  for  steady  state,  self-field  MPD  thrusters[1].  This 
theory  is  capable  of  explaining  both  the  observed  oscillations  as  well  as  the  increased 
erosion  observed  at  Onset.  This  work  represents  the  first  time  that  the  plasma 
discharge  and  electrode  processes  have  been  coupled. 


Erosion  processes  depend  on  a  complex  coupling  between  plasma  discharge 
characteristics,  plasma-wall  interactions,  and  electrode  phenomena.  In  particular, 
erosion  rates  vary  depending  on  whether  the  current  conduction  is  through  localized 
molten  spots,  or  via  a  diffuse  (distributed  mode).  In  either  case,  the  electron  emission 
from  the  electrodes  in  MPD  thrusters  can  be  due  to  field  emission  (quasi-steady 
operation)  or  due  to  thermionic  emission  (steady  state  operation).  It  is  important  to 
emphasize  therefore  that  the  conditions  addressed  in  this  research  are  those 
corresponding  to  steady  state  conditions.  This  means  that  the  theory  and  results 
presented  herein  apply  more  to  the  high  power  steady  state  experiments  being 
conducted  at  the  University  of  Stuttgart,  rather  than  the  quasi-steady  experiments 
being  conducted  at  Princeton  University. 


During  the  period  of  this  work,  a  method  has  been  developed  for  the  systematic 
approximate  determination  of  cathode  temperatures  and  hence,  evaporative  erosion 
rates  in  steady  state,  self-field  MPD  thrusters.  This  theory  uses  the  simplest  level  of 
coupling  between  the  hot  electrodes  and  the  flowing  plasma.  Given  the  propellant 
mass  flow  rate,  total  current,  and  the  geometry,  all  the  relevant  quantities  are  predicted 
for  the  flowing  plasma  as  well  as  electrode  temperatures.  Once  electrode 
temperatures  are  determined,  evaporative  erosion  rates  can  be  determined  from  vapor 
pressure  data[2]: 


log,0  m"  =  7.5  - 


40,500 

T 


(D 


where  T  is  in  degrees  Kelvin,  and  m"is  in  g/cm2/s.  This  model  has  been  compared 

with  steady  state  experiments  at  the  University  of  Stuttgart  on  the  ZT-1  and  DT-2 
thrusters,  as  well  as  with  the  quasi-steady  experiments  on  the  straight  co-axial  and 
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half  scale  flared  anode  thrusters  (HSFAT)  at  Princeton  University.  It  is  found  that  the 
model  correctly  predicts  thrust  and  electrical  characteristics  in  the  quasi-steady 
experiments,  as  well  as  the  observed  cathode  damage  in  the  steady  state 
experiments. 

This  report  in  conjunction  with  three  earlier  annual  progress  reports, 
summarizes  progress  made  under  grant  AFOSR-87-0360. 

2.0  RESEARCH  OBJECTIVES 

The  overall  objective  of  this  work  was  to  develop  the  design  methods  and  tools 
necessary  to  predict  the  limits  of  stable  steady  state  operation  of  the  electrodes  in  MPD 
thrusters.  Given  the  mass  flow  rate,  total  current,  and  geometry,  the  designer  must 
know  not  only  the  thrust  and  electrical  characteristics,  but  also  whether  the  Onset 
condition  has  been  reached  and  whether  the  electrodes  are  operating  under  steady, 
non-molten  conditions.  To  accomplish  this  goal,  simple  and  approximate  models  had 
to  be  found.  These  include  models  of  quasi  one-dimensional  ionizing  MPD  flow, 
electrode-adjacent  boundary  layer  flow,  the  electrode-adjacent  sheath,  and  the 
electrodes  themselves. 

3.0  STATUS  OF  RESEARCH  EFFORT 

The  discussion  below  summarizes  the  highlights  of  the  research  efforts  under 
this  grant.  A  more  detailed  discussion  can  be  found  in  the  appendices  and  earlier 
references!  1,3, 4], 
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Steady  state  operation  nf  MPD  thrusters  has  been  limited  by  the  Onset 
phenomenon.  Onset  is  a  term  that  is  used  to  denote  collectively,  increased  erosion  of 
thruster  components  and  terminal  voltage  oscillations  that  are  seen  to  occur  at  a 
critical  current  for  a  given  mass  flow  rate  and  geometry.  Several  theories  have  been 
proposed  that  can  predict  the  Onset  limit[1 ,3,4-7].  We  focus  here  on  the  back-EMF 
theory  of  Onset[1 ,3]. 


The  idea  of  a  high  back-EMF  being  responsible  for  Onset  phenomena  was  first 
discussed  within  the  context  of  a  one-dimensional,  steady  state,  frozen,  fully  ionized 
flow.  This  model  revealed  that  MPD  flow  was  parametrized  by  the  Magnetic  Force 

Number  S*=B*2/n0Fa*,  where  the  superscript  *  refers  to  quantities  evaluated  at  the 
magnetogasdynamic  sonic  point.  S*  was  also  related  to  the  Onset  parameter  Jc2/m, 


utilized  by  experimentalists  to  describe  Onset.  In  this  theory,  S*  was  found  to  have  an 


upper  limit  for  supersonic  flow  in  the  MPD  thruster,  thus  yielding  a  quantitative  limit  for 
the  Onset  parameter  derived  from  first  principles  by  theory  to  be: 


m 


*  8.52 


W  a 
H  h0k2 


(2) 


where  Jc  is  the  critical  current  at  Onset,  m  is  the  total  propellant  mass  flow  rate,  W  is 
the  effective  channel  width,  H  is  the  effective  channel  height,  a*  is  the  frozen  speed  of 

sound  at  the  sonic  point,  n0  is  the  permeability  of  free  space,  and  k=B/B*  is  the  ratio  of 

magnetic  induction  at  the  inlet  to  the  magnetic  induction  at  the  sonic  point.  Equation 
(2)  was  found  to  give  excellent  agrrem«nt  with  the  experimental  results  of  Malliaris  et. 
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al.[8].  Incorporation  of  finite  rate  ionization  and  recombination  did  not  alter  these 
results  qualitatively,  but  an  explicit  closed  form  expression  for  the  Onset  limit  could  not 
be  found[3]. 

The  back-EMF  theory  has  also  been  compared  with  measurements  on  the 
straight  co-axial  thruster[3],  and  the  half-scale  flared  anode  thruster  (HSFAT)[9]  at 
Princeton  for  quasi-steady  operation.  We  have  shown  previously  that  the  quasi  one¬ 
dimensional  model  predicts  the  electrical  characteristics  for  the  straight  co-axial 
geometry  fairly  well.  Fig.  1  shows  the  thrust  versus  total  current  for  the  HSFAT 
operating  under  quasi-steady  conditions,  it  can  be  seen  that  the  quasi  one¬ 
dimensional  model  compares  extremely  well  with  experiemental  measurements.  In 
fact,  the  agreement  is  even  better  than  the  two-dimensional,  fully  ionized,  frozen,  axi- 
symmetric  model  of  ref[10].  The  transition  from  electrothermal  to  electromagnetic 
acceleration  is  evident  in  the  slight  change  in  slope  visible  in  the  figure  at  a  current 
level  near  14  kA.  The  2-D  model  does  not  appear  to  reproduce  this  fact. 

Although  this  back-EMF  theory  quite  satisfactorily  predicts  the  limits  of  steady 
state  operation,  electrode  phenomena  had  to  be  included  to  study  the  erosion 
processes  in  select  regions  of  the  thruster.  This  is  because  electrode  processes 
together  with  the  plasma-electrode  interactions  influence  erosion.  Further,  no  existing 
theory  was  able  to  explain  the  observed  increase  in  erosion  at  Onset.  In  order  to 
resolve  this  issue,  the  aforementioned  quasi  one-dimensional  model  was  coupled  via 
a  two-temperature  boundary  layer  theory  to  the  electrode-adjacent  sheath.  Analyses 
of  the  sheaths  via  simple  models  revealed  some  new  phenomena: 
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•  Stable,  diffuse  mode  behavior  transitions  into  the  spot  mode  at  a  critical 
value  of  the  sheath  voltage  drop. 

•  This  diffuse  to  spot  transition  is  caused  by  a  thermal  runaway  due  to 
excessive  electron  bombardment.. 

•  The  thermal  runaway  occurs  in  the  low  current  density  regions  of 

the  cathode,  and  the  high  current  density  regions  of  the 
anode. 

The  increased  erosion  at  Onset  can  thus  be  explained  as  follows.  As  the  total  current 
is  increased  for  a  given  mass  flow  rate  and  geometry,  the  back-EMF  climbs  from  a 
small  value  at  the  inlet,  reaches  a  maximum  somewhere  in  the  middle  of  the  channel, 
and  then  decreases  toward  the  exit.  This  causes  the  net  current  density  to  reach  a 
minimum  in  the  middle  of  the  channel,  while  being  large  in  the  exit  region.  From 
current  conservation,  it  can  be  shown  that  the  sheath  voltage  drop  follows  the  same 
trend  as  the  net  current  density[1 1].  Thus,  the  theory  predicts  that  the  middle  portion  of 
the  cathode  (i.e.  away  from  the  inlet  and  exit  regions)  is  prone  to  damage  because  of 
thermal  runaway.  Simultaneously,  when  the  sheath  voltage  drop  reaches  values 
comparable  with  the  energy  of  the  first  excited  state  of  the  propellant  atoms,  the 
electron  distribution  function  is  modified  near  the  electrode,  due  to  two  principal 
groups  of  electrons.  The  first  are  “low”  energy  plasma  electrons,  and  the  second  are 
beam  electrons  emitted  by  the  electrode  and  accelerated  through  the  sheath. 
Interaction  between  these  two  groups  of  electrons  can  lead  to  longitudinal  voltage 
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oscillations  arising  from  the  well  known  beam  instability[12]. 


Experiments  confirm  these  theoretical  predictions.  Barnett  reported  that  for 
uniform  flow  to  the  benchmark  thruster  operating  quasi-steady,  at  Onset,  voltage 
oscillations  were  observed  in  the  exit  regions  of  the  flow[13].  Recent  experiments  at 
the  University  of  Stuttgart  on  high  power  steady  state  MPD  thrusters  showed  severe 
damage  in  the  middle  of  the  cathode  consistently[14],  A  schematic  of  the  damage  is 
shown  for  these  DT-2  and  ZT-1  thrusters  in  Fig.  2  and  Fig.  3.  Both  cathodes  appear  to 
have  exploded  from  within,  at  6500  A  and  8000  A  respectively.  The  back-EMF  theory 
not  only  correctly  predicted  these  current  limits,  but  also  the  temperature  profiles 
shown  in  Fig.  4.  As  expected,  the  centerline  temperature  far  exceeds  the  cathode 
surface  temperature.  Further,  at  Onset,  the  interior  temperature  exceeds  the  melting 
temperature.  When  this  occurs,  pressure  builds  up  inside  the  cathode  driving  crack 
propagation  and  promoting  fissures.  Although  the  back-EMF  theory  correctly  predicts 
this  Onset  and  damage  limit,  it  cannot  describe  the  events  after  the  damage  and 
melting  have  begun  to  occur. 

The  theory  that  has  been  developed  and  discussed  extensively  in  the 
appendices,  enable  the  designer  to  use  quick  tools  for  performance  studies  and 
design  evaluations.  Further  work  is  however  needed  to  accurately  determine  the  level 
of  ionization  in  the  electrode  pre-sheath  regions,  which  can  affect  the  accuracy  of  the 
predicted  erosion  rates  using  the  present  theory. 
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Fig.  1 :  Comparisons  of  total  thrust  versus  total  current  between  simulations  and 

experiment  are  shown  here  for  Princeton's  Half-Scale  Flared  Anode  Thruster  (HSFAT). 


Fig.  2:  Schematic  of  damage  sustained  by  the  Stuttgart  DT-2  thruster  after 

steady  state  operation  up  to  6.5  KA. 


Fig.  3:  Schematic  of  damage  sustained  by  the  Stuttgart  ZT-1  thruster  after 

steady  state  operation  up  to  8  KA. 


Length  along  cathode  (cm) 


Fig.  4:  Shown  here  are  calculated  surface  and  centerline  temperature  profiles 

versus  distance  along  the  cathode  for  the  Stuttgart  ZT*1  thruster.  The  horizontal  line  at 
3680  K  denotes  the  melting  temperature  of  tungsten.  Temperatures  above  this  value 
cannot  be  taken  seriously  from  the  theoretical  predictions,  since  the  theory  does  not 
account  for  phenomena  once  melting  has  occurred.  However,  as  the  model  suggests 
and  as  shown  in  this  figure,  the  centerline  temperature  has  exceeded  the  surface 
temperature  in  the  select  region  thereby  explaining  the  exploded  cathode  observed  in 
the  Stuttgart  test. 
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Appendix  A 


Quasi  One-Dimensional  Ionizing  MPD  Flow 


G.  Lefever-Button  and  V.  V.  Subramaniam 
Department  of  Mechanical  Engineering 
The  Ohio  State  University 
Columbus,  Ohio  43210 
U.S.A. 


Abstract 

Quasi  one-dimensional  models  can  serve  as  useful  testbeds  for  testing  various 
numerical  algorithms  for  integrating  the  types  of  stiff  differential  equations  that 
result  from  including  finite  rate  kinetics  (i.e.  ionization  and  recombination).  They  are 
an  important  prelude  to  more  complicated  multi-dimensional  simulations.  In  this 
paper,  an  earlier  quasi  one-dimensional  model  is  extended  to  study  the  effects  of 
varying  area  on  ionizing  MPD  flow.  Three  popular  numerical  integrators  of  chemical 
rate  equations  (Runge-Kutta,  Gear,  and  Bulirsch-Stoer)  are  compared  for  speed  and 
accuracy.  Based  on  the  tests  reported  here,  the  Bulirsch-Stoer  algorithm  appears  far 
superior.  Results  are  also  presented  for  MPD  channel  flow  with  area  variation  and 
finite  rate  kinetics.  Comparisons  of  the  thrust  computed  using  the  present  quasi  one¬ 
dimensional  model  with  the  recent  experimentally  measured  values  for  the  HSFAT 
(Half-Scale  Flared  Anode  Thruster)  are  excellent.  An  important  outcome  of  this  work 
is  the  finding  that  channel  divergence  has  significant  effects  not  only  on  the  thrust, 
but  also  on  ionization  and  current  density  profiles. 


Nomenclature 


a  Frozen  speed  of  sound 

A  Area 

B  Magnetic  induction 

e  Electronic  charge 

E  Electric  field 

h  Enthalpy  per  unit  mass 

j  Current  density 

J  Total  current 

k  Boltzmann’s  constant 

kf  Ionization  rate  constant 

kr  Recombination  rate  constant 

l  Computed  channel  length 

L  Given  channel  length 

mA  Atomic  mass 


1 


m  Mass  flow  rate 

M  Mach  number 

P  Pressure 

T  Temperature 

u  Velocity 

x  Coordinate  along  channel 

a  Ionization  fraction 

e4  Ionization  potential  of  gas 

Permeability  of  free  space 
p  Mass  density 

a  Electrical  conductivity 

*  superscript  for  quantities  evaluated  at  the  sonic  point 
e  subscript  denoting  exit  values 

i  superscript  denoting  ith  iteration 

i  subscript  denoting  inlet  values 


I.  Introduction 

Earlier  theoretical  works  have  considered  one-dimensional  and  quasi  one-dimensional 
models  of  MPD  flows  (see  for  instance  refs.[l]-[8]).  The  main  motivation  for  such 
drastic  simplifications  has  been  to  study  the  influence  of  individual  phenomena  on 
MPD  flow,  either  as  a  prelude  to  studies  of  complicated  multi-dimensional  plasma 
flows  or  to  arrive  at  a  basic  understanding  of  operating  limits  and  Onset  phenomena. 
However,  there  exists  another  important  use  for  these  models.  Quasi  one¬ 
dimensional  models  that  include  finite  rate  kinetics  can  act  as  test  beds  for  numerical 
algorithms  that  may  be  used  in  multi-dimensional  simulations.  Alternatively,  simple 
models  are  often  very  effective  design  tools.  It  is  with  these  goals  in  mind  that  we 
extend  a  previous  one-dimensional  flow  model  with  finite  rate  kinetics[5],  to  include 
varying-area  channels. 

Finite  rate  kinetics  introduce  numerical  stiffness  that  often  overwhelms  the  speed 
and  even  stability  of  the  best  algorithms[9].  Further,  quasi  1-D  models,  unlike 
multi-dimensional  simulations,  are  complicated  by  the  existence  of  singular  points’ 
such  as  saddle  points,  nodes,  and  spirals  (foci)[10J.  In  this  paper,  the  back-EMF 
theory[4,5]  is  extended  to  include  the  effects  of  area  variation.  Also,  the  best 
presently  available  numerical  techniques  for  handling  stiff  systems  are  compared. 

This  paper  is  organized  as  follows.  The  model  assumptions  and  governing  equations 


are  discussed  next,  followed  by  a  discussion  of  the  numerical  procedures  used  to 
integrate  these  equations.  The  results  for  select  cases  are  then  discussed,  followed 
finally  by  the  summary  and  conclusions. 


II.  Governing  Equations 

We  consider  here  the  equations  of  quasi  1-D,  steady  flow  with  the  electric  field 
perpendicular  to  the  flow  direction.  The  self-generated  magnetic  induction  is  taken 
to  be  perpendicular  to  both  this  electric  field  and  the  flow  direction.  This  is 
tantamount  to  the  assumption  of  negligible  Hall  effect.  Further,  we  assume  the  MPD 
plasma  to  be  composed  of  electrons,  neutrals,  and  singly  charged  ions.  Finite  rate 
kinetics  are  included  in  the  overall  collisional  reactions  of  electron-impact  ionization 
and  three-body  recombination: 


+A  ~e~  +A*  +  e~ 


(1) 


The  ionization  rate  kf  and  recombination  rate  kr  for  argon  are[5]: 

48  x  10~9exp(-e(/kT) 

k,= - - -  (cm  /s) 

f  T\ 5.556  xlO11)3'2 

*,  =  4xl<r9r9/2  (cm6 /s)  * 


where  e(  is  the  ionization  potential  of  the  propellant.  The  conservation  equations  are: 
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energy: 
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Gauss’  law: 


Ampere’s  law: 
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Although  viscous  and  two-temperature  effects  are  known  to  be  important  11],  they 
are  neglected  here.  The  existence  of  two  distinct  flow  regimes,  one  dominated  by 
ohmic  heating  and  the  other  by  the  electromagnetic  body  force,  leads  to  the 
magnetogasdynamic  sonic  or  choking  condition  [5].  This  choking  condition 
determines  the  electric  field  at  the  sonic  point: 

£*  =-amB‘  +{p,  +  p2  +  p,)1/2  <“> 
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where 


p,  =  —  a*V 

1  16 


,  and 

2  omA‘dx 
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„  p*fl*  *i  da 

and  *  denotes  quantities  evaluated  at  the  sonic  point,  i.e.  where  the  flow  speed  equals 
the  local  frozen  speed  of  scund[5].  This  sonic  or  choking  condition  has  been  the 
subject  of  much  confusion  [12,13],  but  a  rigorous  treatment  including  proofs  can  be 
found  in  ref.[10]. 


III.  Numerical  Procedure 

Before  proceeding  with  a  discussion  of  the  solution  procedure  of  the  system  (3) 
through  (10),  we  briefly  describe  three  integrators  of  systems  of  ordinary  differential 
equations.  These  integrators  are  (1)  4th  order  Runge-Kutta  with  variable  step  size, 
(2)  Gear’s  predictor-corrector  method[14],  and  (3)  Bulirsch-Stoer  technique  with 
Richardson  extrapolation[15]. 

In  the  4th  order  Runge-Kutta  method,  functional  values  at  the  initial  point,  two 
intermediate  points,  and  a  trial  endpoint  are  used.  Derivatives  are  then  evaluated 
at  each  of  these  four  points,  and  a  weighted  average  of  these  values  is  used  to 
calculate  the  end  point  value. 

The  Bulirsch-Stoer  method  uses  rational  function  extrapolation  (Richardson 
extrapolation)  in  order  to  extrapolate  interim  solutions  to  that  corresponding  to  a  zero 
step  size.  This  is  achieved  by  successively  bisecting  the  original  interval  and 
subsequent  subintervals,  and  integrating  using  the  midpoint  method  (2nd  order 
Runge-Kutta)[16].  In  contrast  to  these  two,  Gear’s  method  uses  many  initial  and/or 
previous  points  to  predict  the  functional  value  at  the  next  point.  A  corrector  equation 
specially  formulated  for  stiff  equations  is  then  used  to  iterate  on  the  functional  value 
at  the  next  point. 

Although  all  three  numerical  methods  are  capable  of  handling  stiff  equations,  they 
vary  considerably  in  the  handling  of  the  sonic  singularity.  Consequently,  these 
techniques  were  first  tested  on  the  quasi  1-D  isentropic  flow  equations  of  classical  gas 
dynamics.  This  test  isentropic  flow  problem  consisted  of  an  area  variation  given  by' 
A(x)=1.98x2-3.96x+2  (m2),  channel  length  of  2  meters,  upstream  stagnation  pressure 
of  1.01x10s  N/m2,  upstream  stagnation  temperature  of  273  K,  and  p  =1.7798759259 
Kg/m  .  A  comparison  of  the  three  numerical  methods  in  terms  of  the  number  of  right 
hand  side  evaluations,  number  of  bisections  to  evaluate  the  correct  inlet  velocity,  and 
total  elapsed  CPU  time  (run  on  a  VAX  8500),  is  shown  in  Fig.  1.  In  Fig.  2,  the 


overall  accuracy  of  the  three  methods  is  compared.  It  can  be  seen  that  the  best 
compromise  between  speed,  efficiency,  and  accuracy,  is  the  Bulirsch-Stoer  algorithm. 
This  is  evident  when  all  three  methods  are  compared  for  speed  while  holding  the  total 
error  constant  to  within  the  same  order  of  magnitude.  In  this  case,  the  Bulirsch- 
Stoer  algorithm  is  seen  to  be  most  efficient,  followed  by  the  Gear  method,  and  then 
by  the  4th  order  Runge-Kutta  method  with  variable  step  size. 


The  system  (3)  through  (10)  may  be  simplified  since  the  density,  electric  field,  and 
enthalpy  may  be  computed  directly  from: 


pio4  =  constant 

(12) 

EA  =  constant 

(13) 

2  p  mA 

(14) 

This  leaves  the  following  four  equations  to  be  integrated: 

ax 

(15) 

da  _k/PaO-ay  *rP2°3 
*  mAu  m\u 

(16) 

5  Pu  dA  2  *i  da  5  u  n.  2-. 
du  3  pA  dx  3  mA  dx  3  p  3p 

dx  u2  -  a2 

(17) 

—  =  Ba(E-  uB)  -  p«— 
dx  dx 

(18) 

The  above  system,  along  with  the  associated  boundary  conditions,  constitutes  a  two- 
point  boundary  value  problem.  At  the  nozzle  inlet,  the  mass  flow  rate,  magnetic 
induction,  and  ionization  fraction  are  specified.  The  nozzle  area  variation  as  a 
function  of  position  x  is  given,  as  is  the  total  thruster  length.  The  velocity  and 

(q 


electric  field  at  the  inlet  must  be  solved  for,  such  that  the  downstream  flow  exiting 
the  thruster  is  supersonic  and  the  computed  length  matches  the  specified  thruster 
length.  The  temperature  at  the  inlet  is  constrained  by  the  value  of  the  electric  field, 
which  is  itself  bound  by  the  sonic  point  choking  condition.  However,  the  location  of 
the  sonic  point  and  hence  the  value  of  the  electric  field  at  that  point  are  not  known 
a  priori.  The  correct  inlet  values  can  only  be  determined  after  the  sonic  point  is 
successfully  crossed.  The  algorithm  developed  thus  entails  integrating  through  the 
channel  repeatedly,  each  time  refining  the  inlet  conditions  using  bisections.  Since  the 
sonic  point  is  a  mathematical  singularity,  a  "jump"  algorithm  is  needed  to  cross  it. 
This  is  discussed  next. 

An  iterative  procedure  is  used  to  cross  the  sonic  point.  This  proceeds  as  follows.  Let 
the  adjacent  point  on  the  subsonic  side  of  the  sonic  point  be  denoted  the  subscript  1, 
and  the  adjacent  point  on  the  supersonic  side  by  the  subscript  2.  Let  M2=2-M1.  The 
temperature  at  the  pre-sonic  point  Tx  is  used  as  an  initial  approximation  for  T2. 
Then  u2  is  computed  iteratively  from: 

=  M2  +  ajf-1))  '  <19> 

where  the  superscript  i  refers  to  the  ith  iteration.  The  jump  distance  is  then: 

the  new  temperature  is  computed  from: 

r2(0  =  T.  +  —  (21) 

2  1  dx  x,Xl 

and  the  new  ionization  fraction  is: 

a<<>  =  o1  +  —  Ax(0  (22) 

2  1  dxx.Xl 

Equations  (19)  through  (22)  are  placed  in  a  loop  where  they  are  repeatedly  evaluated 
until  the  differences  in  temperature  and  ionization  fraction  between  two  successive 
iterations  is  smaller  than  some  specified  tolerance.  In  the  present  calculation  a 
tolerance  of  10  6  was  specified.  Thus  the  Mach  number,  velocity,  temperature,  and 
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ionization  fraction  on  the  supersonic  side  of  the  sonic  point  are  all  calculated,  along 
with  its  coordinate: 

*2°  =  x  j  +  Ax(<)  ^ 

The  other  properties  across  the  sonic  point  are  then  determined  using  the 
conservation  equations. 

The  integration  of  the  system  (15)  through  (18)  is  achieved  starting  from  the  subsonic 
inlet.  Since  not  all  initial  conditions  are  known,  usually  a  manual  shooting  method 
is  used  in  which  guessed  initial  estimates  are  refined  after  integration  past  the  sonic 
point.  Rather  than  relying  on  such  an  ad-hoc  shooting  method,  an  automatic 
algorithm  for  integrating  the  governing  equations  has  been  developed.  To  accomplish 
this,  a  sequence  of  tests  have  been  designed  which  determine  which  inlet  variable  to 
modify  and  in  what  direction  to  modify  it.  Some  of  these  tests  are  performed  after 
a  complete  integration  through  the  channel,  while  others  are  done  in  a  specific  region, 
such  as  within  the  jump  routine.  Still  others  are  performed  after  every  integration 
step.  These  diagnostic  tests  are  summarized  in  Table  I. 

As  mentioned  previously,  the  inlet  temperature  is  constrained  somewhat  by  the 
electric  field.  It  was  found  empirically  that  this  constraint  was  much  weaker  for 
lower  values  of  total  current.  It  was  also  found,  however,  that  the  solution  was  not 
sensitive  to  the  inlet  temperature.  Figure  3  shows  a  comparison  of  the  resulting 
temperature  distributions  for  total  current  J  =  8.276  kA,  and  inlet  temperatures 
•  ranging  5500  <  T;  <  11,500  K  for  the  flared  geometry  of  the  HSFAT  thruster,  which 
will  be  discussed  in  the  following  section.  For  the  case  of  T;  <  5500  K  convergence 
could  not  be  obtained,  and  at  T4  >  11,500  K  the  flow  is  already  supersonic  at  the  inlet. 
It  can  be  seen  that,  beyond  the  sonic  point,  the  results  do  not  change  for  different 
inlet  temperatures.  The  computed  thrust  varied  up  to  a  maximum  of  7%. 

The  numerical  algorithm  presented  here  has  been  applied  to  the  solution  of  the 
governing  equations  for  a  constant  area  channel.  These  results  have  been  compared 
to  existing  solutions  for  constant  area  channels[5].  The  agreement  is  found  to  be  good 
to  at  least  the  third  decimal  place  in  all  the  variables.  In  the  next  section,  results  are 
presented  for  the  simplest  of  geometries,  i.e.  linearly  varying  channel  height. 

IV.  Results 

Before  discussing  the  varying  area  cases,  we  will  briefly  revisit  the  constant  area 
geometry.  The  effect  of  different  mass  flow  rates  on  the  profiles  of  ionization  fraction 
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along  the  length  of  the  channel  is  shown  in  Fig.  4,  for  mass  flow  rates  ranging  from 
0.003  kg/s  to  0.006  kg/s.  The  magnetic  induction  at  the  inlet  (x=0)  is  0.159,  the  cross 
sectional  area  is  0.001  m2,  the  length  is  20  cm.,  and  the  inlet  ionization  fraction  is 
0.05.  It  can  be  seen  from  Fig.  4  that  the  density  has  a  profound  effect  on  the  rate  of 
ionization.  From  equation  (7),  it  can  be  seen  that  the  recombination  term  in  the  rate 
equation  for  a  varies  as  p2  in  contrast  to  the  ionization  term  which  varies  as  p. 
Thus,  for  higher  mass  flow  rates,  and  hence  higher  densities,  a  actually  decreases 
along  the  channel  until  near  the  exit,  where  the  temperature  increases.  This  effect 
of  density  on  the  ionization  fraction  has  not  been  previously  reported[5-8,12]. 

Let  us  now  consider  a  channel  with  an  area  variation  given  by  A(x)=0.001+x(8A), 
where  x  is  in  meters  and  A(x)  is  in  m2.  The  total  channel  length  is  taken  to  be  20 
cm.  so  that  xe[0.0,0.2].  The  four  cases  reported  here  are  5A=0  (constant  area), 
5A=0.0175,  5A=0.04,  and  8A=0.07.  The  fixed  quantities  are  inlet  area  Ai=10  cm2, 
mass  flow  rate=0.006  kg/s,  B^O.159,  <^=0.05,  and  length=20  cm.  The  profiles  of  u, 
T,  j,  a,  B,  and  M  are  shown  in  Figs.  5-10.  Qualitatively,  each  profile  is  similar  to 
those  discussed  previously [5].  However,  the  flared  varying  area  geometry  is  seen  to 
be  very  different  in  several  ways  when  compared  to  the  constant  area  geometry. 
First,  the  most  obvious  effect  of  varying  area  is  apparent  in  the  velocity  distributions. 
The  exit  speeds  and  hence  specific  impulses  are  1.7,  2.0,  and  2.2  times  larger  than 
the  constant  area  case  for  area  ratios  of  4.5,  9,  and  15  respectively.  Second,  the 
ionization  fraction  shows  the  same  behavior  seen  in  Fig.  4  for  a  constant  area  channel 
with  the  higher  mass  flows.  The  flared  geometry  has  a  dramatic  effect  on  the 
ionization  fraction.  As  the  flare  angle  is  increased,  the  value  of  a  is  seen  to  increase. 
This  is  due  to  the  aforementioned  decrease  in  density  for  the  same  mass  flow  rate  but 
different  flare  angles.  Third,  the  magnetic  induction  shows  the  effects  of  increasing 
back-EMF  as  the  flare  angle  is  increased  for  a  fixed  mass  flow  rate  and  total  current. 
Finally,  the  current  density  rise  at  the  exit  usually  seen  in  the  constant  area 
channels[5]  is  diminished  greatly  by  the  flared  geometry. 

We  now  turn  our  attention  to  a  modified  version  of  the  original  King  anode[17]. 
Thrust  measurements  have  recently  been  made  on  the  Half-Scale  Flared  Anode 
Thruster  (HSFAT)[18J.  Comparisons  between  these  HSFAT  measurements  and  two- 
dimensional  axi-symmetric  frozen,  fully  ionized  flow  modelling  results  have  also 
recently  been  performed[19].  The  HSFAT  is  essentially  the  same  shape  as  the' 
original  King  anode[17],  but  much  shorter  in  axial  length  (10  cm.  versus  20  cm.). 
Figures  11  and  12  show  comparisons  of  thrust  and  thrust-to-total  current  squared, 
between  the  present  work,  experimental  measurements  of  ref.[18],  and  the  numerical 
results  of  ref.[19].  As  can  be  seen,  the  agreement  is  excellent.  Furthermore,  there 
is  a  subtle  inflection  point  around  11  kA  in  the  experimental  results  shown  in  Fig. 


11  that  is  correctly  reproduced  in  the  present  results.  This  inflection  point  is  due  to 
the  transition  from  electrothermal  to  electromagnetic  acceleration.  Taking  the 
electrothermal  contribution  to  the  thrust  to  be  proportional  to 

L 

Jj(E-uB)A  dx  (24) 

o 

and  the  electromagnetic  contribution  as 

L 

f[u-(jxB)]Adx 

o 

and  comparing  the  ratio  of  electromagnetic  to  electrothermal  contributions  gives  the 
graph  of  Figure  13. 

In  addition,  there  is  a  dramatic  increase  in  thrust  beyond  approximately  21  kA.  The 
present  quasi  one-dimensional  predictions  end  around  16  kA  because  the  flow  is 
nearly  fully  ionized  over  a  significant  length  inside  the  thruster  (see  Fig.  14).  We 
conclude  therefore  that  the  dramatic  increase  in  thrust  at  currents  above  20  kA  is 
due  to  the  presence  of  doubly  charged  ions,  which  we  have  not  considered  here.  It 
is  also  interesting  to  note  the  typical  current  density  and  back-EMF  profiles  shown 
in  Figs.  15  and  16.  Also  plotted  in  these  figures  for  reference  are  the  corresponding 
curves  for  a  linear  flared  anode  (i.e.  arm  varying  linearly  from  inlet  to  exit).  The 
back-EMF  is  highest  in  the  region  where  the  electric  field  is  highest,  and  lower  where 
the  electric  field  is  lower.  This  is  in  contrast  to  the  straight  coaxial  thruster  where 
the  electric  field  is  constant  and  the  back-EMF  reaches  a  maximum  in  the  region 
approximately  midway  between  inlet  and  exit.  Such  nozzle  contouring  can  therefore 
help  control  Onset  due  to  an  excessive  back-EMF[4,5]. 


V.  Summary  &  Conclusions 

An  efficient  algorithm  utilizing  the  Bulirsch-Stoer  method  for  integration  of  quasi 
one-dimensional  ionizing  MPD  flow  has  been  presented.  Solutions  have  been 
compared  with  previous  solutions  for  constant  area  channels,  and  have  revealed 
features  not  reported  in  previous  studies[5-8,12].  This  work  has  revealed  that  for  a‘ 
constant  channel  area  and  total  current,  the  exit  ionization  fraction  increases  with 
decreasing  mass  flow  rate.  Furthermore,  for  a  fixed  mass  flow  rate  and  total  current, 
substantial  increases  in  ionization  fraction  and  specific  impulse  may  be  achieved  by 
using  a  flared  geometry.  Too  much  of  a  flare  or  abrupt  expansion,  however,  can  lead 
to  Onset  by  the  back-EMF  mechanism[5].  Comparison  of  the  present  model  with 


measurements  on  the  HSFAT  show  excellent  agreement.  This  suggests  that  thruster 
geometries  can  be  optimized  for  specific  operating  points  (i.e.  total  current  and  mass 
flow  rate)  to  maximize  thrust  and  specific  impulse,  using  the  present  model. 
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Figure  and  Table  Captions 


Figure  1:  Comparison  of  integrator  performance  characteristics  for  isentropic  test 
case. 

Figure  2:  Comparison  of  overall  integration  errors  for  isentropic  test  case. 

Figure  3:  Plot  of  temperature  profiles  for  varying  inlet  temperatures  in  the  flared 
channel  of  the  HSFAT  thruster  (see  section  IV  for  details). 

Figure  4:  Ionization  fraction  profiles  for  a  constant  area  channel,  A=0.001  m2, 
Bj=0.159,  L=20  cm.,  at  various  mass  flow  rates. 

Figure  5:  Velocity  profiles  for  exit/inlet  area  ratios  of  1  (constant  area),  4.5,  9,  and 
15,  with  Ai=0.001  m2  and  20  cm.  length. 

Figure  6:  Temperature  profiles  for  linearly-varying  areas,  ApO.OOl  m2  and  20  cm. 
length. 

Figure  7:  Current  density  profiles  for  linearly-varying  areas,  ApO.OOl  m2  and  20 
cm.  length. 

Figure  8:  Ionization  fraction  profiles  for  linearly-varying  areas,  Ai=0.001  m2  and  20 
cm.  length. 

Figure  9:  Magnetic  field  profiles  for  linearly- varying  areas,  A—0.001  m2  and  20  cm. 
length. 

Figure  10:  Mach  number  profiles  for  linearly-varying  areas,  Ai=0.001  m2  and  20  cm. 
length. 

Figure  11:  Comparison  of  experimentally  measured  and  calculated  values  of  thrust 
for  10  cm.  HSFAT  channel. 

Figure  12:  Comparison  of  experimentally  measured  and  calculated  values  of 
thrust/J2  for  10  cm.  HSFAT  channel. 

Figure  13:  Comparison  of  the  ratio  of  electromagnetic/electrothermal  contributions 
for  various  total  currents. 


Figure  14:  Comparison  of  ionization  fraction  profiles  in  10  cm.  HSFAT  channel. 

Figure  15:  Comparison  of  HSFAT  current  density  profile  with  that  of  a  linearly- 
varying  channel  with  the  same  At  and  Ae. 

Figure  16:  Comparison  of  back-EMF  of  HSFAT  with  that  of  a  linearly-varying 
channel  with  the  same  Ai  and  Ac. 

Table  I:  Tests  performed  and  actions  taken  for  automatic  integration  of  MPD 
channels.  Variables  subscripted  1  are  previous  points,  those  subscripted  2  are 
evaluated  at  the  present  location,  and  the  subscript  T  denotes  the  throat,  or 
minimum  area  location. 
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Abstract 

A  great  amount  of  research  has  been  performed  in  self-field 
Magnetoplasmadynamic  (MPD)  thruster  research  during  the  past  decade.  In 
contrast  to  previous  years,  activity  in  the  theoretical  and  modelling  areas 
have  increased.  This  paper  will  review  the  developments  in  the  theory  of 
self-field  MPD  thrusters,  concentrating  mainly  on  one-dimensional  and 
quasi  one-dimensional  theories,  for  which  exact  or  approximate  solutions 
can  be  obtained. 

a  Ionization  fraction 

Cj  Ionization  potential 

Heavy  particle  viscosity 

k  Ratio  of  magnetic  induction  at  inlet  to 

the  magnetic  induction  at  the  chokinq 
point 

M0  permeability  of  free  space 

p  Mass  density 

o  Electrical  conductivity 

i  Subscript  denoting  inlet  quantity 

e  Subscript  denoting  exit  quantily 

Superscript  denoting  quantity  evaluated 
at  the  sonic  or  choking  point 

I.  Introduction 

During  the  past  ten  years,  there  has  been  a 
great  deal  of  activity  in  self-field 
magnetoplasmadynamic  (MPD)  thruster 
research.  This  paper  will  attempt  to 

summarize  the  developments  in  the  theory  of 
these  thrusters  and  compare  theory  to 
available  experiments.  Models  are  now 
capable  of  predicting  l-V  curves  from  first 
principles.  The  conditions  under  which 
electrodes  erode  are  beginning  to  be 
understood.  Quantitative  prediction  of  the 
destructive  ’onset"  phenomenon  also  appears 
to  be  possible  now. 

Early  work  in  the  field  of  MPD  thrusters 
has  been  summarized  by  Jahn[1).  The 
technology  ol  MPD  thrusters  has  been 


Nomenclature 

a  Frozen  speed  of  sound 

B  Magnetic  induction 

Ce  Mean  thermal  speed  of  electrons 

C^  Mean  thermal  speed  of  heavy  particles 

E  Electric  field 

h  Enthalpy  per  unit  mass 

H  Channel  height 

j  Current  density 

J  Total  current 

k  Boltzmann's  constant 

L  Channel  length 

m  Mass  flow  rate 

mA  Atomic  mass 

m0  Electron  mass 

n„  Electron  or  charged  particle  number 

density 
P  Pressure 

QAA  Total  momentum  transfer  cross  section 

for  atom-atom  collisions 
q  A  Total  momentum  transfer  cross  section 

for  ion-atom  collisions 
q  I  Total  momentum  transfer  cross  section 

for  ion-ion  collisions 
T  Temperature 

T0  Electron  temperature 

Th  Thrust 

Tj_(  Heacy  particle  temperature 

u  Axial  component  of  velocity 

W  Channel  width 

x  Coordinate  along  channel 
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recently  reviewed  by  Myers  el  al[2).  This 
review  will  concenlrale  on  developments  in 
the  theory  ol  MPD  thrusters.  As  two 
dimensional  models  are  still  in  their  early 
stages  o(  development,  emphasis  will  be 
given  to  quasi-one-dimensional  models  and 
their  implications  (or  thruster  per.ormance 
and  future  thruster  design. 

Magneloplasmadynamic  (MPD)  thrusters 
use  electric  and  magnetic  fields  to 
accelerate  a  plasma  to  high  speeds,  commonly 
20  to  45  km/s.  This  acceleration  is  achieved 
by  the  action  of  jxB  lorces  in  a  steady  flow. 
The  magnetic  field  may  be  generated  by  the 
current  itself,  "self-field."  Alternatively, 
the  magnetic  induction,  B.  may  be  externally 
applied.  Sell-field  thrusters  generally 
require  instantaneous  input  powers  of  the 
order  of  1  MW.  Applied  field  thrusters  can 
operate  with  powers  as  low  as  100  kW.  Over 
the  past  ten  years,  research  has  concentrated 
on  s el  I- field  thrusters  and  this  review  will 
concentrate  on  this  type. 

A  sample  configuration  for  a  self-field 
thruster  is  shown  in  Figure  1.  The  geometry 
is  axisymmelric.  The  propellant,  typically  a 
rare  gas  and  most  often  argon,  is  injected 
through  holes  in  the  back  plate.  A  current  is 
drawn  between  the  center  cathode  and  outer 
anode.  The  current  induces  an  azimuthal 
magnetic  field.  The  resulting  jxB  force  is 
largely  axial,  accelerating  the  propellant  out 
ol  the  Ihrusler. 

Self-field  thrusters  tend  to  have  higher 
efficiencies  as  the  current  driving  them 
increases.  Specific  impulse  also  increases 
with  current.  This  increase  in  performance 
has  a  well-know  limit:  above  a  certain 
current  level,  instabilities  and  erosion  set  in. 
This  is  collectively  called  "onset"  and  the 
current  at  which  this  happens  is  called  the 
critical  current  or  Onset  current.  Much 
theoretical  ellort  has  been  devoted  towards 
explaining  this  phenomenon,  and  hence  will 
also  be  discussed  in  this  review. 

This  paper  is  organized  as  follows. 
Early  Mow  models  are  discussed  in  the  nexl 
section,  followed  by  a  discussion  of  theories 
ol  Onset  Electrode  and  boundary  layer 
models  are  then  discussed,  tollowed  by  the 
current  status  ol  analyses  of  electrode 
behavior  under  steady-stale  MPD  operation. 


It.  Flow  Models 

White  the  flow  in  an  MPD  thruster  may  be 
complex,  there  is  one  simple  feature:  the 
magnetic  body  force  is  the  strongest  force 
operating  on  the  fluid.  From  this  starting 
point,  useful  information  can  be  derived.  The 
magnitude  of  the  magnetic  body  force  depends 
only  on  the  current  distribution  and  not 
directly  on  the  fluid  properties.  By  assuming 
current  profiles  and  integrating  the  body 
force  over  the  flow  volume,  the  thrust,  Th. 

can  be  derived.  If  the  current  lines  are  radial 
between  a  cylindrical  cathode  of  radius  rc 

and  an  anode  of  radius  rg,  as  shown  in  Figure 
2a,  the  resulting  thrust  is: 


where  J  is  the  total  current  and  p0  is  the 

permeability  of  free  space.  If  the  current  is 
carried  by  the  conical  section  ol  the  cathode 
as  shown  in  Figure  2b,  the  thrust  is 
determined  from: 


A  remarkable  feature  of  these  equations  is 
that  they  do  not  depend  on  the  path  that  the 
current  follows  between  the  electrodes(1  ]. 

Notice  that  the  geometry  in  Figure  2b 
produces  a  larger  thrust  for  the  same  current 
J  than  the  geometry  of  Figure  2a.  This  does 
not,  however,  mean  that  one  geometry  is  more 
efficient  than  the  other.  This  level  of  theory 
cannot  tell  us  the  voltage  required  to  drive 
the  current  J  or  hence  the  electrical  input 
power  used  to  generate  the  thrust  Th. 

The  next  level  of  theory  considers  one¬ 
dimensional  (or  quasi  one-dimensional)  flow 
through  the  thruster.  A  thruster  geometry 
which  can  be  analyzed  as  one-dimensional  is 
shown  in  Figure  3.  This  plane-parallel 
geometry  avoids  radial  variations  in  lields 
that  occur  when  rg»rc.  The  behavior  is 

otherwise  analogous.  For  example,  by 
conservation  of  momentum,  the  thrust  for 
this  geometry  can  be  found  to  be: 

t  "oj2  H 

Th  ”  2  W  ^ 

This  equation  is  very  similar  in  lorm  to  those 
above  In  (act,  equation  (3)  can  by  derived 
from  equation  (1)  by  identifying  FUrg-  rc  and 

W=2ar c  and  considering  1-1  small.  Since  the 


2 


IIhusI  T(l  is  Ihe  product  of  the  mass  (low 

rale,  m,  and  Ihe  exhaust  speed, u0,  the  above 
equation  also  determines  exhaust  speed: 

^  H 


For  the  case  o(  an  applied  (ield  thruster, 
the  gas  dynamics  ol  the  channel  (low  were 
first  studied  in  1958  by  Sears  and 
Resler[3.4j.  Their  study  mapped  Ihe  various 
regions  ol  operation  including  such  important 
effects  as  MPD  choking.  Further  discussion  of 
this  theory  can  be  found  textbooks{5,6). 

In  the  1970's,  attention  shifted  to  self¬ 
field  thrusters.  Self-field  thrusters  offered 
engineering  simplicity  by  not  requiring  large 
magnets.  Inspired  by  a  Ph.O.  thesis  by 
Martinache[7).  King[8.9],  developed  a  model  of 
self-field  thrusters  that  could  successfully 
compare  to  experiment  In  King's  model,  the 
thruster  is  modeled  as  a  one-dimensional 
flow  subjected  to  ohmic  heating  and  the 
magnetic  body  force  (See  Table  I). 

To  understand  the  flow,  it  is  necessary 
to  know  how  Ihe  various  phenomena  affect 
the  flow's  Mach  number.  This  is  summarized 
in  Table  II  and  is  not  always  intuitive(IO).  At 
Ihe  inlet  ol  a  MPD  thruster,  ohmic  healing  is 
dominant  and  this  drives  Ihe  How  towards 
Ihe  sonic  point  (M=1).  As  the  flow  speed 
increases,  the  magnetic  body  force  becomes 
increasingly  important.  Under  exactly  the 
right  condition,  ohmic  healing  lifts  the  Mach 
number  to  1  and  the  magnetic  body  force  then 
drives  the  Mach  number  still  higher.  The 
condition  under  which  this  happens  is  called 
the  choking  condition,  and  a  rigorous, 
mathematical  discussion  of  this  condition 
can  be  found  in  Ref. (11).  For  an  ideal 
monatomic  gas,  Ihe  choking  condition  is: 

E  =  |  a4  B'  (5) 

where  a  is  the  speed  of  sound,  E  is  the 
eleclric  field,  B  is  Ihe  magnetic  induction  and 
Ihe  superscript  *  indicates  that  these  are 
evaluated  at  the  choking  point.  This  condition 
must  be  satisfied  or  (he  flow  cannot 
accelerate  from  subsonic  to  supersonic 
speeds.  This  generally  happens  far  upstream 
in  the  thruster  as  shown  in  Figure  1. 

One  might  naively  assume  that,  il  the  flow 
is  introduced  into  the  thruster 
supersonically,  then  the  choking  condition 
would  be  irrelevant  This  is  generally  not  the 
case  Immediately  after  the  How  is  injected, 
it  is  still  subjected  to  strong  ohmic  heating 
From  Table  II,  ohmic  healing  lowers  the  Mach 


number  of  supersonic  fiow.  Since  Ihe 
magnitude  of  the  ohmic  heating  generally 
greatly  exceeds  the  stagnation  enthalpy  of 
the  injected  propellant,  the  How  is  quickly 
driven  to  Mach  1.  In  order  to  accelerate  again 
to  supersonic  speeds,  the  choking  condition 
must  still  be  met.  The  case  of  supersonic 
propellant  injection  was  extensively  studied 
by  King[9J. 

King  recognized  that  this  model  had 
important  implications  (or  thruster  design. 
This  led  to  the  development  of  the  King 
anode[8j. 

To  clarify  the  relationship  between  this 
theory  and  earlier  work,  Ihe  equations  { 1  )- 
(3).  (4)  for  thrust  can  be  derived  from  the 
momentum  equation  given  in  Table  I  by 
neglecting  pressure. 

There  is  another  restriction  on  the 
electric  field.  From  consideration  of  Ohm's 
law  under  usual  MPD  thruster  conditions,  it  is 
found  that  the  back-EMF  must  never  exceed 
the  applied,  electric  field,  E: 

uBsE  (6) 

This  restriction  does  not  replace  the  choking 
condition.  Ralher,  it  must  be  obeyed  in 
addition  to  the  choking  condition,  equation 
(5).  This  presents  an  interesting  and 
important  issue.  Let  us  combine  Ihe 
inequality  (6)  with  the  choking  condition: 

uB  S  | a  *  B  *  (7) 

Somewhere  in  the  middle  region  of  the 

thruster,  the  product  uB  reaches  a  maximum. 
This  region  is  observed  experimentally  as  a 
region  of  low  current  density.  Let  us  call  the 
values  of  u  and  B  at  this  point  um  and  Bm, 

respectively.  Now,  rewrite  inequality  (7)  as: 

5  0* 

'  23’  B^  (8) 

The  magnetic  induction  at  Ihe  choking  point, 
B\  is  within  10%  of  the  magnetic  induction 
at  the  thruster  inlet.  The  magnetic  induction, 
Bm,  where  the  back-EMF  is  a  maximum  is 

usually  0  6  to  0.7  limes  the  value  at  the  inlet. 
Thus,  their  ratio  can  be  approximated  as  1.5. 
Substituting  in  this  value: 

um  <  3  75a‘  (9) 

The  right-hand  side  above  is  fairly  constant, 
a’  is  determined  by  the  temperature  at  the 
choking  point  and  plasma  temperatures  are 

known  to  vary  little.  The  left-hand  side 
varies  strongly  with  current.  From  equation 

(4),  the  exhaust  speed  scales  as  J^/  m.  Thus. 

equation  (9)  appears  to  place  a  limit  on  J2/  m 
Let  us  estimate  um  as  2/3  of  u0.  Then 


3 


combining  equalion  (4) 
gives: 


2  /'0J2  N 

3  2m  W 


3.75a  * 


with  equalion  (9) 


PO) 


Rearranging  this  gives: 

5 1 1 25n7"  c  * 

m  H  /'o 

This  equalion  expresses  the  restriction  on 

o  * 

J  / m  that  comes  from  combining  the  choking 
condition,  equation  (5),  with  the  limit  on 
back-EMF,  uBsE.  The  above  derivations 
included  several  approximations.  II  the 
equations  on  Table  I  are  solved  exactly,  the 
result  is: 


~  S  8  52 


W  a‘ 
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where  k  is  the  ratio  of  the  magnetic  induction 
at  the  inlet  to  the  magnetic  induction  at  the 
choking  point.  This  typically  has  a  value  of 
about  1.1(9).  This  explanation  of  onset  first 
appeared  in  Rel(12). 


field  below  the  peak  values  that  the 
equilibrium  theory  would  allow.  For  a 
monatomic  gas  with  finite  rate  ionization, 
the  choking  condition  is: 


2  I  mA  dx 


(14) 


where  o  is  the  ionization  traction,  tj  is  the 
ionization  potential,  m^  is  the  atomic  mass, 

and  j‘  is  the  current  density  evaluated  at  the 
choking  point.  It  is  apparent  that  the 
ionization  rate  at  the  choking  point  raises  the 
electric  held  above  the  value  it  would  have 
for  an  ideal  monatomic  gas.  This  allows  the 

thruster  to  operate  at  higher  J2/m  than  the 
ideal  gas  theory  would  permit.  How  much 
higher  depends  on  the  ionization  traction  and 
temperature  at  the  choking  point.  To 
determine  this  currently  requires  a  full 
numerical  model.  With  such  a  model,  it  has 
been  possible  to  quantitatively  predict  onset 
currents,  as  shown  in  Figure  5  and  Figure 
6[13). 


It  should  be  mentioned  that  the  rote  of 
inlet  flow  conditions  in  altering  the 
occurrence  of  back-EMF  onset  is  contained 
entirely  in  the  k  factor.  This  quantity  is 
restricted  to  a  small  range  about  1.1.  Inlet 
conditions  are  extensively  discussed  in  Ref. 
(9J 

The  above  discussion  was  based  on  the 
ideal  gas  choking  condition,  equalion  (5).  For 
quantitative  results,  non-ideal  effects  must 
be  considered,  principally  ionization.  If  one 
considers  equilibrium  Mow,  the  choking 
condition  becomes: 

f3> 

The  quantity  in  square  brackets  is  5/2  lor  an 
ideal  monatomic  gas.  For  equilibrium  argon, 
its  value  is  shown  in  Fig  4.  It  is  seen  that 
the  equilibrium  theory  and  ideal  monatomic 
gas  theory  agree  only  at  low  temperatures, 
<5000K.  For  higher  temperatures,  the 
equilibrium  theory  predicts  much  higher 
electric  Melds  These  higher  electric  fields 
would  have  the  ellect  ol  delaying  back-EMF 

onset  to  larger  values  ol  J2/m  The 

magnitude  of  the  increase  in  film  can  be 
estimated  by  replacing  5/2  in  equation  (7) 
with  the  values  given  in  Fig.  4  and  then 
continuing  the  derivation  as  before. 


For  simplicity,  most  ol  the  models  above 
have  assumed  a  one-temperature  plasma.  The 
plasmas  in  MPD  thrusters  are  more  accurately 
described  by  two-temperatures.  A  fairly 
complete  simulation  of  a  two-temperature 
MPD  thruster  with  ionizing  argon  has  been 
performed  by  Richley{14).  The  simulation 
was  capable  ol  modeling  fully  unsteady  How. 
A  sample  temperature  profile  is  shown  in 
Figure  7. 

Shoji  and  Kimura  have  developed  a  theory 
ol  two-temperature  (low  in  MPD 
thruslers(15J.  They  also  extended  the  model 
to  include  dissociating  molecular  gases  and 
performed  calculations  on  hydrogen.  As  (or 
ionization,  nonequilibrium  dissociation  alters 
the  choking  condition,  raising  the  allowed 
electric  held.  They  restricted  their  study  to 
the  region  between  choking  and  exit  regions 
and  parametrically  varied  the  conditions  at 
choking.  For  fixed  thermodynamic  conditions 
at  the  choking  point  and  for  fixed  total 
current,  they  varied  the  mass  (low.  They 
numerically  found  a  minimum  mass  flow.  i.e. 

maximum  fit  m  ,  for  which  the  thruster  had  a 
finite  length.  This  i:  exactly  the  behavior 
expected  from  the  above  discussion  ol  back- 
EMF  Onset  and  choking  condilions.  This  is 
the  first  calculation  ol  this  effect  in  a  Iwo- 
lemperature  model. 


Unfortunately,  the  equilibrium  assumption 
is  not  accurate  Real  gases  have  finite 
ionization  rales  and  this  reduces  the  electric 


Models  of  self-lield  MPD  thrusters  have 
greatly  advanced  over  the  past  ten  years. 
Beginning  with  simple  estimates  for  thrust. 
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models  are  now  able  to  predict  detailed 
temperature  and  speed  distributions.  The 
models  described  in  this  section  have, 
however,  been  restricted  to  quasi  one¬ 
dimensional  Hows.  This  limits,  lor  example, 
their  ability  to  describe  How  in  the  Princeton 
"benchmark’  thrusters  which  have  an  anode 
lip  protruding  into  the  (low.  They  have, 
however,  quantitatively  described  Hows  in 
straight  channel  coaxial  thrusters  as  well  as 
thrusters  with  the  diverging  King  anode{8,9, 
13).  More  recent  comparisons  between  this 
quasi  one-dimensional  model  and 
experimental  results  at  Princeton  on  the 
Hall-Scale  Flared  Anode  Thruster  (HSFAT) 
show  excellent  agreement  as  can  be  seen 
Irom  Fig.  8[16J.  This  suggests  that  quasi  on- 
dimensional  models  may  serve  as  very  useful 
tools  for  design  and  optimization. 

III.  More  Theories  o(  Onset 

The  phenomenon  ol  "onset"  has  attracted  a 
great  deal  of  experimental  and  theoretical 
interest  Several  theoretical  explanations  ol 
onset  have  been  advanced  and  more  than  one 
have  experimental  evidence  to  support  them. 
It  may  be  possible  that  onset  does  not 
represent  a  single  phenomenon  but  several 
phenomena  which  all  occur  at  high  current 
levels.  The  back-EMF  theory  ol  onset  was 
discussed  above  because  it  is  an  integral  part 
ol  quasi-one-dimensional  MPD  llow  theory. 
Some  of  the  other  theories  ol  onset  will  be 
discussed  here. 

Several  theories  locus  on  plasma 
conditions  near  the  anode.  Specifically, 
consider  a  plasma  with  electron  density  n0 

and  mean  thermal  speed  Ce=(OkT0/  m^)1 
The  random  thermal  current  density  is 
ithermareneCe/4  ,f  ,,ie  ne*  curren|  density 
between  the  plasma  and  the  anode  is  (ess 
than  jlherma(,  an  anode  sheath  forms  with  a 

potential  that  discourages  the  llow  of 
electrons  to  the  anode.  This  potential  is  a 
negative  voltage  drop  and  is  observed  most  of 
the  time  during  low  current  operation  ol  MPO 
thrusters  For  high  current  operation,  several 
effects  combine  to  change  this  condition. 
First,  lor  high  current  operation  the  current 
density  is  simply  higher.  Secondly,  at  high 
current,  the  llow  is  accelerated  to  higher 
speeds  and  hence  lower  densities.  Thirdly,  af 
higher  currents,  the  Hall  effect  becomes 
stronger  and  this  tends  to  drive  the  plasma 
Irom  the  anode  towards  the  thruster 
centerline  These  latter  two  effects  tend  to 


reduce  the  random  thermal  current  density. 
•  thermal'  At  sulficiently  high  current,  the 
required  net  current  density  may  exceed  the 
random  thermal  current.  Under  this 
condition,  the  anode  sheath  voltage  drop  is 
expected  to  change  signs  and  become  positive 
in  order  to  help  draw  more  current.  This 
effect  is  sometimes  called  anode  mass 
starvation. 

This  change  in  sign  of  the  anode  sheath 
voltage  drop  was  analyzed  by  Baksht, 
Moizhes,  and  Rybakov  in  1974  using  a  simple 
quasi-one-dimensional  model[17j.  They  took 
this  change  in  sign  as  their  onset  criterion. 
Several  refinements  to  this  model  have  been 
made  by  Heimerdinger[18],  Korsun  considered 
a  similar  mass  starvation  elfect  in  a  two 
dimensional  modeH19}.  Hugel  measured  the 
anode  sheath  voltage  drop  in  a  thruster  and 
found  that  onset  was  associated  with  the 
change  in  sign[20). 

The  reason  why  a  change  in  sign  ol  the 
anode  sheath  should  cause  onset  was  not 
clear  because  many  discharges  operate  in  a 
stable  steady  state  with  anode  sheaths  ol 
both  signs.  Two  hypotheses  exist  to  explain 
this.  Kuriki  and  Onishi  suggest  that  when  the 
anode  tall  is  positive,  ions  accelerated  in  the 
shealh  may  produce  spultering{21  J. 
Alternatively,  Shubin  found  that  anode  mass 
starvation  is  associated  not  just  with  anode 
sheath  reversal  but  also  with  plasma 
instabili!ies[22].  He  suggested  that  these 
instabilities  may  cause  onset. 

II  onset  were  solely  due  to  anode  mass 
starvation,  the  experimental  remedy  would  be 
clear:  inject  more  mass  near  the  anode.  This 
in  fact  was  done  as  part  ol  a  very  detailed 
set  of  onset  experiments  performed  by 
Barnett[23).  He  found  that,  under  some 
conditions,  onset  was  associated  with 
instabilities  near  the  anode.  However,  he 
also  found  that  injecting  mass  near  the  anode 
did  not  prevent  the  occurrence  ol  onset 
Rather,  onset  still  occurred  but  was 
associaled  wilh  instabilities  in  olher  regions 
of  the  thruster.  He  found  that,  by  controlling 
mass  injeclion,  he  could  cause  the  region  ol 
instability  to  move  around  the  thruster  Irom 
anode  to  cathode.  This  appears  to  indicate 
that,  while  anode  mass  starvation  may  be  one 
onset  mechanism,  it  is  certainly  not  the  only 
onset  mechanism. 
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one  moves  radially  toward  the  walls  and 
away  from  the  core  flow,  the  gas  becomes 
less  ionized.  This  is  because  of  wall-driven 
recombination  reducing  the  number  density  of 
charged  particles.  In  these  wall-adjacent 
regions,  the  electron  temperature  tends  to 
remain  fairly  constant,  while  the  heavy 
particle  (ions  and  neutrals)  temperature 
decreases  from  the  core  flow  value  down  to 
very  nearly  the  wall  or  electrode 
temperature.  The  decrease  In  ionization 
fraction  as  one  approaches  the  walls  causes 
the  viscosity  to  be  determined  by  the  ion- 
neutral  total  momentum  transfer  cross 
section  (which  is  typically  several  orders  of 
magnitude  smaller  than  Qj  j.  Thus,  the 

viscosity  changes  by  orders  of  magnitude 
going  from  the  core  to  the  walls.  If  the 
ionization  fraction  in  the  near-wall  regions 
drops  much  below  -0.05,  then  dominates 

the  viscosity.  Since  MPD  flows  are  less  than 
fully  ionized  over  much  of  the  channel  length, 
they  are  also  viscous. 


V.  Steady-State  MPD  Operation 

It  is  important  to  understand  erosion  and 
lifetime  issues  in  MPD  thrusters.  Because  of 
the  high  power  required,  most  laboratory 
experiments  on  MPD  thrusters  have  involved 
operation  for  only  a  few  milliseconds.  In 
space,  by  contrast,  it  is  likely  that  MPD 
thrusters  will  be  operated  in  steady-state 
The  erosion  and  lifetime  issues  are  very 
different  between  these  two  modes. 


IV.  Electrode  and  Boundary  Layer 
Models 

A  key  issue  in  MPD  thruster  research  has 
been  the  electrode  erosion  rales  and 
mechanisms.  Between  the  bulk  plasma  flow 
and  the  wall,  there  is  a  fluid  mechanical 
boundary  layer  and  an  eleclroslatic  boundary 
layer  (sheath).  These  layers  govern  the 
plasma-electrode  interactions. 

Between  the  bulk  plasma  How  and  the 
eleclrodes.  there  are  (hermal.  viscous,  and 
electrostatic  boundary  layers.  Modeling  of 
these  layers  is  very  important  to  both 
thruster  performance  and  life.  They  affect 
performance  because  they  determine  the 
viscous  drag.  They  affect  life  because  they 
determine  plasma  conditions  adjacent  to  Ihe 
electrodes  and  hence  influence  electrode 
operating  temperatures  and  heat  fluxes. 

fn  the  boundary  layer,  the  flow  speed 
changes  from  the  free  stream  value  to  a 
small  value  just  above  the  electrode  surface. 
Several  other  plasma  properties  also  change, 
including  ionization  fraction,  heavy  particle 
temperature,  and  plasma  density.  Not  all 
these  properties  change  monotonically. 
Compressible  boundary  layer  theory  for 
single-temperature  fluids  has  been 
summarized  by  Schlichling[24j. 

I|  is  only  recently  that  a  quantitative 
model  of  the  MPD  boundary  layer  has  become 
available(25j  This  boundary  layer  analysis, 
Ihe  first  such  lor  an  MPD  thruster  showed  the 
important  effects  of  ionization  fraclion  on 
viscosity,  and  hence  on  boundary  layer 
growth  The  following  expression  lor 
viscosity  can  be  derived  from  kinetic 
theory[25j: 

nH  -  mAcf|aq(*  (,-ajqA  *(i-BtX3M*BqA)  n<3) 
where  C is  Ihe  mean  thermal  speed  of  Ihe 

heavy  particles,  <r  is  Ihe  ionization  fraction, 
C(  ^  is  (he  total  momentum  transfer  cross 

section  (or  ion-neulral  collisions,  Q^a  is  (lie 

total  momentum  transfer  cross  section  for 
neutral-neutral  collisions,  and  Qj  j  is  the 

total  momentum  transfer  cross  section  for 
ion  ion  collisions  In  Ihe  core  flow,  i.e  away 
from  Ihe  walls,  the  plasma  is  nearly  fully 
ionized  so  that  tr  is  near  1.  In  (his  region,  Ihe 
momentum  transfer  is  determined  primarily 
the  ion-ion  collisions,  as  can  be  seen  from 
equation  (16)  by  selling  o-1  How  over,  as 


When  operated  for  a  few  milliseconds,  Ihe 
plasma  flow  has  lime  to  reach  steady-state. 
The  electrodes,  however,  do  not  have  time  to 
reach  thermal  steady-state.  As  a 
consequence,  the  erosion  mechanism  in  pulsed 
("quasi-steady”)  operation  and  steady-slate 
operation  are  very  different. 

For  Ihe  cold  electrodes  of  quasi-steady 
operation,  il  Is  experimentally  observed  that, 
at  high  current  levels,  visible  hoi  spots  form 
on  Ihe  electrode  surfaces.  These  hot  spots 
can  be  explained  as  thermal  instabilities  of 
Ihe  plasma-elecirode  interaction.  Theories 
have  recently  been  developed  to  describe 
these  e(fecls[26-28J  The  nature  of  these 
instabilities  differs  between  the  cathode  and 
anode. 
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For  Ihe  long  life  demanded  of  space¬ 
worthy  equipment,  the  erosion  rates  of  spot¬ 
mode  current  conduction,  >1(ig/C,  are 
generally  unacceptable.  On  hot  cathodes, 
current  can  be  emitted  (hermionically.  This 
enables  erosion  rates  to  drop  down  to  the 
order  of  1  ng/C(29J. 

An  important  issue  is  the  stability  of 
diffuse  (thermionically)  emitting  cathodes. 
Two  recent  papers  have  analyzed  thermal 
stability  limits  of  diffuse  mode  operation  of 
cathodes  and  anodes(26,27).  The  limits  have 
been  found  to  be  very  different  for  the  two. 
Both  stability  limits  involve  thermionic 
emission  and  sheath  theory. 

On  anodes,  it  is  found  that  areas  of  low 
current  density  are  generally  stable.  In 
regions  of  high  current  density,  a  runaway 
heating  of  the  electrode  can  occur.  The 
runaway  is  driven  by  thermionic  emission. 
While  thermionic  emission  is  undesirable  on 
an  anode,  it  happens  unavoidably  at  the  high 
temperatures  (eg.  3000K)  of  steady-slate 
operation.  While  thermionic  emission  cools 

Ihe  anode,  it  has  another  indirect  effect.  The 
more  thermionic  emission  that  there  is,  the 
more  plasma  current  that  must  be  drawn  to 
the  electrode.  The  net  eflect  is  heating  since 
the  plasma  electrons  are  much  hotter  (e  g. 
10.000K)  than  the  anode.  The  result  of 
heating  Ihe  anode  is  Still  more  thermionic 
emission  and  then  still  more  heating.  The 
anode  spots  observed  in  Ihe  high  current 
regions  of  steady-stale  thrusters  may  be 
evidence  for  this  mechanism. 

On  cathodes,  thermionic  emission  is 
essential  for  diffuse  mode  operation.  This 
emission  must  be  adequate  to  conduct  the 
eleclricly  in  the  high  current  density  regions, 
such  as  the  lip  of  Ihe  cathode.  Since  the  mid- 
section  of  the  cathode  is  usually  nearly  as 
hot  as  the  lip  (or  possibly  more  so),  it  also 
has  thermionic  emission.  Its  emission, 
however,  becomes  excessive  as  onset  is 
approached.  Near  onset,  the  current  density 
in  Ihe  middle  of  the  thruster  drops  rapidly 
due  to  an  increasing  back-EMF.  The  result  is 
that  plasma  electrons  must  enter  the  cathode 
to  replace  the  electrons  emitted 
thermionically  This  again  leads  to  a  thermal 
runaway  effect 

The  recent  theories  that  have  been 
developed  for  cathodes  and  anodes  operating 
at  steady  state  predict  several  effects  that 
may  be  observable  in  constant  cross  sectional 


area  geometries.  First,  at  high  total 
currents,  current  densities  are  highest  near 
the  inlet  and  exit  regions  while  being 
smallest  in  the  middle  region  due  to  high 
back-EMF.  Second,  thermal  runaway 
destabilizes  Ihe  diffuse  mode  operation  at 
high  current  densities  at  the  anode.  The  same 
mechanism  destabilizes  the  diffuse  mode  at 
the  cathode  at  low  current  densities.  This 
means  that  portions  of  Ihe  anode  near  the 
exit,  and  portions  of  the  cathode  near  the 
middle  are  particularly  susceptible  to 
melting.  Recent  steady  slate  experiments  on 
straight  coaxial  thrusters  at  Stuttgart  appear 
to  corroborate  this  prediction[28]. 

VI.  Summary 

This  review  article  has  discussed 
recent  developments  in  Ihe  theoretical 
understanding  of  self-field  MPD  thrusters. 
The  focus  has  been  primarily  on  one¬ 
dimensional.  quasi  one-dimensional,  and 
partially  two-dimensional  models.  Ongoing 
efforts  in  two  and  three  dimensional 
simulations  of  steady  MPD  flow[30-32]  are 
likely  to  extend  our  present  understanding  of 
MPD  How  dynamics  and  behavior.  Electrode 
processes  are  also  begining  to  be  understood, 
and  when  fully  coupled  with  flow  models, 
should  provide  valuable  insight  in  the  design 
of  experiments  and  interpretation  of 
experimental  results. 
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Mass:  pu  *  F  »  constant 


Momentum: 


F 


«  constant 


Energy: 


—  »  constant 


State:  h  •  h(P,p) 

Ohm's  Law:  j  »  o(E-uD) 

*  /  ,  dO 

Ampere  s  Law:  —  *  \ij 

dx 

Tabic  I.  Equations  of  motion  for  a 
self-Held  onc-tcmperatorc  quasi-onc- 
dimensional  Mt’D  flow. 


n.  2:  Idealized  ntodcls  with  (a)  uniform 

.Jli.t  r-tirrrnl  — *  i",'‘  "liliral 


Fig-  3:  Onc-rlintcnsioihil  thruster 

gcoinetry. 


Phenomena 

Effect  on 

Mach  number 
in  subsonic 

flow 

F.ffcct  on 

Mach  number 
in  supersonic 
flow 

Ohmic  hentinz 

increases 

decreases 

Accelerating 

body  force 

decreases 

increases 

Friction 

increases 

decreases 

Wall  cooline 

decreases 

increases 

Endothcmiic 

reactions 

decreases 

increases 

. 

Table  If:  Various  plienouicna  cause  lltc  flow  Macli  numticr 
to  change  toward  or  away  front  I. 


Fig.  4:  The  electric  field  determined  by 
the  mngnetogasdynamic  choking  condition, 
is  plotted  here  against  the  temperature  at 
the  choking  point,  for  the  two  cases  of 
equilibrium  and  frozen  flow  (from  ref.(  1 2 J). 
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Fig.  5:  Electric  field  versos  total  current, 
for  a  total  inass  (Tow  of  3  g/s.  The 
uppermost  point  on  |he  solid  curve  is  lire 
theoretical  Onset  limit  (from  ref  f)3|). 


Fig.  t>:  Electric  field  versus  total  current, 
for  a  total  mass  now  of  f>  g/s.  The 
uppermost  point  on  the  solid  curve  is  the 
theoretical  Onset  limit  (from  ref.(l3l). 


Fig.  7:  Heavy  particle  and  electron 
temperatures  in  a  20  cm.  long  MPD 
thruster,  from  rcf.fMJ. 


Fig.  R:  Comparison  of  thrust  versus  total 
current  between  quasi  one-dimenswna) 
model,  two-dimensional  axisymmclric 
model,  and  experimental  measurements  for 
Ihc  IISFAT  lltmsicr  (from  ref-t  1 6]). 
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Abstract 

A  method  is  outlined  for  the  determination  of  cathode  temperatures  and 
hence,  evaporative  erosion  rates  using  the  simplest  level  of  coupling 
between  the  electrode  and  flowing  MPD  plasma.  Given  the  mass  flow  rate, 
total  current,  and  geometry,  a  previously  developed  quasi  one¬ 
dimensional  model  is  used  to  generate  profiles  of  all  relevant  quantities 
in  the  axial  direction.  These  are  then  used  as  boundary  conditions  for  a 
boundary  layer  model,  which  in  turn  generates  the  necessary  boundary 
conditions  for  the  electrode-adjacent  sheath  and  electrode  thermal 
analysis.  The  steady  state  heat  conduction  equation  is  then  solved 
numerically  to  estimate  cathode  temperatures.  The  coupling  between  the 
flowing  plasma  and  the  electrode  is  modelled  via  the  boundary  conditions 
so  that  solutions  are  tractable.  A  simple  version  of  this  technique  is  then 
used  to  compare  model  predictions  with  experimental  observations  of 
cathodes  in  steady  state  MPD  thrusters  operating  at  high  current  levels. 
The  model  correctly  predicts  the  observed  damage  to  the  ZT-1  and  DT-2 
thrusters,  and  shows  the  inter-relations  between  back-EMF  onset,  erosion, 
and  terminal  voltage  oscillations. 


a  Speed  of  sound 

B  Magnetic  induction 

C  Mean  thermal  speed 

e  electronic  charge 

H  Channel  height 

j  Net  current  density 
J  Total  current 

k  Boltzmann's  constant 

kc  Thermal  conductivity  of  cathode 
material 

L  Thruster  length 

mA  Atomic  mass 

q  Heat  flux 

qend  Heat  ,lux  at  cathode  tip 

%aseHeat  flux  at  base 
qsurface  Heat  flux  at  the  cathode 
surface 

S  Magnetic  force  number 

T  Temperature 

u  Velocity 

W  Channel  width 

x  Coordinate  along  cathode  length 
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a  Ionization  fraction 

q  Viscosity  of  heavy  particles 

H  Permeability  of  free  space 

p  Mass  density 

o  Electrical  Conductivity 

*  Superscript  indicating  evaluation 
at  the  sonic  point 

i  Subscript  denoting  inlet  value 

I.  Introduction: 

The  behavior  of  self-field  MPD 
thrusters  at  high  values  of  total  current 
for  a  given  propellant  mass  flow  rate 
and  thruster  geometry,  has  been  a 
subject  of  considerable  interest. 
Particularly,  the  phenomenon 
encountered  at  high  currents  known  as 
Onset  has  received  much  attention  due 
to  its  significant  influence  on  the 
performance  of  these  devices.  Onset  is 
a  term  used  to  denote  collectively, 
increased  erosion  of  thruster 
components  and  terminal  voltage 
oscillations  that  are  seen  to  occur 
during  operation  at  high  current 
levels[1,2).  The  focus  of  this  paper  is 
to  show  that  the  conditions  leading  to 
voltage  oscillations  and  electrode 
erosion  are  indeed  linked,  and  to  provide 


» 


a  basis  for  quantitative  prediction  of 
erosion  rates  at  and  below  the  Onset 
limit. 

Several  theories  have  been  successful 
at  predicting  an  Onset  limit[3-6] . 
Schrade  et.  al.  have  attempted  further 
to  explain  the  increased  erosion 
observed  at  Onset[4].  They  have 
analyzed  the  behavior  and  stability 
limits  of  a  single  current  carrying 
channel,  and  have  predicted  the  critical 
current  at  which  such  a  channel 
becomes  unstable  subsequently  forming 
other  neighboring  spots  (molten, 
concentrated  regions  of  high  current 
density  on  an  electrode).  Although  the 
erosion  rate  from  a  single  spot  can  be 
predicted  with  such  a  model, 
generalization  to  an  overall  erosion 
rate  from  an  electrode  surface  is  an 
insurmountable  task. 

Recently,  it  has  been  shown  that 
cathode  and  anode  behavior  are  strongly 
controlled  by  a  combination  of 
variables,  notably  the  local  current  and 
number  densities[7,8).  A  thermal 
runaway  mechanism  triggered  by 
bombardment  of  energetic  plasma 
electrons  shown  to  operate  at  the 
cathode  in  regions  where  the  local 
current  density  was  low  and  the  local 
number  density  was  high(8].  The  same 
mechanism  was  shown  to  be  operative 
at  the  anode  as  well,  but  in  regions  of 
high  current  density  and  low  number 
density(7].  In  this  paper,  it  is  shown 
that  these  conditions  are  consistent 
with  the  Back-EMF  theory  of  Onset,  and 
that  this  theory  is  capable  of  explaining 
oscillations  as  well  as  increased 
erosion  at  Onset.  Also,  a  methodology 
is  outlined  whereby  evaporative  erosion 
rates  can  be  predicted  purely  as  a 
function  of  the  global  parameters  such 
as  total  mass  flow  rate,  total  thruster 
current,  and  thruster  geometry.  Such  an 
approach  is  deemed  a  useful  tool  to  the 
designer. 

This  paper  is  organized  as  follows.  The 
Back-EMF  Onset  theory  is  reviewed 
first  in  the  following  section,  followed 
by  a  detailed  discussion  of  a  simple 
model  that  couples  the  MPO  plasma  flow 
and  electrode  processes.  The  numerical 


results  from  this  model  are  presented, 
followed  by  a  discussion  of  the 
relationship  between  Back-EMF  Onset 
and  erosion,  and  Back-EMF  Onset  and 
oscillations,  respectively.  Finally,  the 
findings  of  the  present  study  are 
summarized  and  conclusions  drawn. 


LL _ Back-EMF  Onset  Revisited: 

The  idea  of  a  high  Back-EMF  being 
responsible  for  Onset  phenomena  was 
first  discussed  within  the  context  of  a 
one-dimensional,  steady  state,  frozen, 
fully  ionized  flow[6].  Although 
restrictive  in  its  assumptions,  this 
simple  model  revealed  several 
important  facts.  First,  MPD  flow  was 
parametrized  by  the  Magnetic  Force 

Number  S  -B  2/p0 Fa  ,  where  the 

superscript  *  denotes  quantities 
evaluated  at  the  magnetogasdynamic 
sonic  or  choking  point.  Second,  there 
was  shown  to  be  an  upper  limit  on  S* 
for  sustaining  supersonic  flow  in  the 
thruster.  Finally,  S*  was  shown  to  be 

related  to  the  Onset  parameter  Jc2/m, 


utilized  by  experimentalists  to  describe 
Onset.  A  limit  on  S'  thus  translated  to 

a  limit  on  J  2/m,  and  this  limit  was 


derived  to  be: 


m 


s  8.52 


W  a 
H  /V2 


0) 


Equation  (1)  showed  excellent 
agreement  with  the  data  of  Malliaris  et. 
al.  on  quasi-steady  MPD  thrusters[1]. 

The  measured  dependence  of  Jc2/m  on 

the  propellant  atomic  mass  as  well  as 
on  the  thruster  geometry  are  correctly 
predicted  by  (1)[6).  Subsequently,  the 
assumption  of  frozen,  fully  ionized  flow 
was  removed[9J.  This  meant  that  closed 
form  analytical  solutions  could  not  be 
obtained  for  MPD  channel  flow. 
However,  numerical  solutions  of  the 
governing  equations  showed  an  upper 
limit  on  S'  analogous  to  that  given  by 
equation  (1).  This  upper  limit  was 

interpreted  to  be  the  Onset  limit,  and 
showed  agreement  with  the 
experimentally  measured  Onset  limits 
in  quasi-steady,  straight  coaxial 
thrusters[10].  Physically,  Back-EMF 
Onset  is  caused  by  a  conflict  between 
the  electric  field  necessary  to  draw  all 
the  applied  current  (as  per  Ohm's  law ) 
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and  that  amount  of  the  electric  field 
necessary  to  sustain  supersonic  flow  in 
the  thruster  (as  dictated  by  the 
magnetogasdynamic  sonic  or  choking 
condition).  Both  constraints  on  E  could 
be  met  at  low  currents,  but  beyond  a 
critical  value  of  the  current  (or  S*), 
both  requirements  could  not  be 
simultaneously  satisfied  at  steady 
state.  For  a  more  recent  discussion  of 
this  theory,  see  ref .[1 1  ]. 

Although  this  Back-EMF  theory  can  quite 
satisfactorily  predict  the  limits  of 
steady  operation  of  self-field 
thrusters,  electrode  phenomena  must  be 
included  to  study  erosion  processes  in 
select  regions  of  the  thruster.  This  is 
because  electrode  processes  together 
with  the  plasma-electrode  interactions 
influence  erosion.  It  is  to  these  that 
we  now  focus  our  attention. 

UL _ A  Simple  Model: 

The  structure  of  MPD  flow  has  been 
examined  previously  within  the  context 
of  a  one-temperature  core  flow,  with  a 
two-temperature  boundary  layer, 
including  finite-rate  ionization  and 
recombination[12).  These  models  have 
revealed  that  the  MPD  flow  is  highly 
viscous  with  entry  lengths  on  the  order 
of  a  few  centimeters[12].  The  primary 
reason  for  this  high  viscosity  is  due  to 
the  effects  of  decreasing  ionization 
fraction  on  the  viscosity[12]: 

As  one  proceeds  from  the  flow 
centerline  toward  the  electrodes 
(walls),  the  ionization  fraction  and 
charged  particle  number  density  drop  due 
to  wall-driven  recombination.  The 
consequences  of  this  can  be  seen  from 
equation  (2)  when  a  varies  from  a 
number  near  1  towards  a  value  near  0. 
The  dominant  cross  section  for 
momentum  transfer  changes  from  that 
due  to  Coulomb  interactions  in  the  core 
flow,  to  that  due  to  ion-neutral 
collisions  in  the  boundary  layer.  Since 
these  two  cross  sections  typically  differ 
by  orders  of  magnitude,  the  viscosity 
gains  increasing  importance  in  the  near¬ 
wall  regions  of  the  flow[12).  It  must  be 
pointed  out  that  several  authors  have 
earlier[  13.14]  and  more  recently(15] 
included  two-temperature  effects,  but 
primarily  in  the  axial  direction.  This  has 
an  effect  on  the  boundary  layer  but  is 
relatively  weak  when  compared  to  the 
effect  of  changing  a,  since  tin  varies  as 


the  square  root  of  T^  and  thus  the 

boundary  layer  thickness  then  only 

-r  1  /4 
varies  as  TH'  . 


With  electrode  phenomena  in  mind,  let 
us  now  consider  the  following  model. 
Under  the  assumption  of  axisymmetry, 
the  steady  state  temperature 
distribution  in  the  cathode  satisfies  the 
heat  conduction  equation,  which  in  polar 


coordinates  is: 


a2  i2 

^k-T)  *  « 


=  0 


(3) 


where  j  is  the  current  density  through 
the  cathode,  o  is  the  electrical 
conductivity,  and  kc  is  the  thermal 

conductivity.  The  boundary  conditions 
are: 


-! 


=  0 


-k, 


-I 

'd  r[-r 


-k 


-k 


H' 

cdx 

*T| 


k=0 


<?xl<=L 


"^surface 

"  %ase  or 
^end 


T(x  =0,  r) 


=Tt 


base 


(4) 

(5) 

(6) 
(7) 


where  qsurface  is  taken  to  be  positive 


when  flowing  into  the  surface  at  r*rc, 
%ase  is  P°si,ive  when  flowing  aul  of 
the  surface  at  x=0,  and  qend  is  positive 

when  flowing  aul  of  the  surface  at  x-L. 
The  condition  (4)  is  a  symmetry 
condition  valid  only  for  a  solid  cathode. 
Since  the  centerline,  r*0,  is  a  line  of 
symmetry,  it  is  an  adiabat.  With  these 
considerations,  the  maximum 
temperature  should  be  expected  to 
occur  somewhere  along  the 
centerline,  i.e  along  r=0.  Because  of 
the  Neumann  boundary  conditions  (4) 
through  (7),  there  is  a  constraint  that 
must  be  met  in  order  for  a  steady  state 
to  exist: 


*[0  ,  r_  dx  +  J jlc  - —  rdrdx 

JMsurface  c  1  0  <7 

0  0 

-  f  %ase  rdr  *  f  %nd rdr  (8) 

0  0 

Equation  (8)  is  an  overall  energy 
balance  on  the  cathode  and  states  that 
the  heat  removed  from  the  cathode 
base,  qbase,  must  sa,isfV  '*  in  order  to 
maintain  a  steady  state. 
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The  surface  boundary  at  r-rc  couples 

the  solid  cathode  to  the  adjacent 
plasma.  The  surface  energy  flux 
consists  of  ion  and  electron 
bombardment,  surface  electron 
emission,  and  radiative  energy  exchange 
with  the  anode.  This  surface  and  sheath 
model  has  been  discussed  in  a  previous 
paper(8j.  qsurfacs  is  given  by: 

Surface  (r“rc  x’T)  “  iW  — “J 

-*•¥) 

‘  1 1 14  ‘  "’anode)  ,9) 

where  the  cathode  sheath  voltage  drop 
Vc  is  determined  from  overall  current 

conservation^].  The  plasma  electron 
current  density  je  is  also  dependent  on 

V_  and  can  be  a  dominant  source  of 
c 

heating  under  MPD  conditions[7,8].  As 
Vc  decreases,  je  increases 

exponentially  leading  to  increased 
surface  heating.  The  cathode  tip  is 
assumed  to  radiate  to  cold  space,  so 
that: 

Wr>  -  ecrSBT4(r'x-L)  po) 

A  method  can  now  be  outlined  to 
determine  cathode  surface 
temperatures.  A  similar  approach  can 
be  used  to  determine  anode 
temperatures  as  well.  As  in  any  design 
or  experimental  situation,  the  total 
current,  propellant  mass  flow  rate,  and 
geometry  are  assumed  to  be  specified. 
Given  these  quantities,  a  suitable  flow 
model  may  be  used  to  provide  axial 
profiles  of  the  plasma  temperature, 
ionization  fraction,  current  density,  and 
velocity  along  the  length  of  the 
thruster[9,16].  Next,  a  two- 
temperature  boundary  layer  model  can 
be  used  to  predict  charged  particle 
number  densities,  current  densities, 
electron  and  heavy  particle 
temperatures  outside  the  sheath 
edge[1 2]  An  appropriate  sheath  model 
can  then  be  utilized  to  calculate  the  ion 
number  densities,  electron  number 
densitites,  sheath  voltage  drop, 
electron  temperatures,  and  heavy 
particle  temperatures  along  the 
electrode  surface[8].  These  provide  the 
necessary  information  to  determine  the 
right  hand  side  of  equation  (9).  The 
system  (3)  through  (10)  is  then  solved 
numerically  to  calculate  the 


temperature  profiles  within  the 
electrodes  as  well  as  on  the  surface. 
Knowing  the  electrode  surface 
temperature,  evaporative  erosion  rates 
can  be  determined  from  vapor  pressure 
data.  In  this  manner,  plasma  discharge 
and  electrode  processes  are  coupled.  It 
must  be  pointed  out  at  the  outset,  that 
the  method  just  outlined  is  only  and 
approximate  solution  procedurefor 

considering  the  fully  coupled  plasma 
flow  and  electrode  processes.  The  fully 
coupled  problem  is  extremely  trdious  to 
solve,  even  numerically. 

In  the  following  section,  the 
approximate  method  just  outlined  is 
applied  to  study  the  ZT-1  steady  state 
thruster  at  the  University  of 
Stuttgart[17], 

IV.  Numerical  Results 

Recently,  there  have  been  reports  of 
tests  conducted  at  the  University  of 
Stuttgart  on  steady-state,  straight, 
coaxial,  self-field  MPD  thrustersjl  7], 
The  results  of  these  tests  showed 
severe  cathode  damage  in  the  middle 
as  opposed  to  the  tip  (see  Figures  1  and 
2).  It  appears  that  the  cathodes  melted 
and  exploded  from  within  suggesting 
that  the  internal  temperatures  far 
exceeded  the  surface  temperatures  at 
this  location.  The  authors  of  ref.  [17] 
report  that  they  believe  the  damage  to 
be  caused  by  ohmic  heating  in  selected 
regions  of  the  thoriated  tungsten 
cathode.  However,  this  phenomenon  can 
and  will  occur  at  Onset  according  to  our 
theory,  and  is  intimately  linked  to  the 
recently  reported  thermal  runaway 
mechanism[7,8]. 

The  method  discussed  in  the  previous 
section  can  be  applied  to  the  ZT-1  and 
DT-2  thrusters  used  in  the  Stuttgart 
experiments.  Since  specific  details 
were  not  reported  for  the  DT-2 
thruster,  we  focus  on  the  ZT-1  thruster. 

This  thruster  had  a  cathode 
approximately  18  cm.  long  and  1.8  cm 
diameter.  The  thruster  consisted  of  a 
straight  constant  area  channel  15  cm. 
long,  with  three  segmented  anode 
sections  accounting  for  the  latter  9  cm. 
of  the  channel  length.  The  argon 
propellant  was  introduced  at  an  axial 
location  of  15  cm.  before  the  exit,  at  a 
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rate  of  2g/s.  This  thruster  according  to 
ref. [17],  was  operated  at  steady  state 
up  to  about  8000  A.  This  information 
regarding  the  geometry,  propellant 
mass  flow  rate,  and  total  current, 
provides  enough  information  to  obtain 
electrode  temperature  distributions 
using  the  method  of  section  III.  An 
additional  simplification  is  made  here. 
The  boundary  layer  analysis  is  omitted 
so  that  the  core  flow  quasi  one¬ 
dimensional  axial  profiles  are  directly 
applied  as  boundary  conditions  to  the 
electrode  thermal  analysis. 

Given  the  total  current,  channel 
geometry,  and  propellant  mass  flow 
rate,  the  governing  equations  of  quasi 
one-dimensional  MPD  flow  including 
ionization  and  recombination,  are 
solved  using  a  method  previously 
reported[9].  This  analysis  yields 
profiles  of  all  the  relevant  quantities, 
figures  3  -  8  show  the  profiles  of  B,  u. 
T,  a,  the  back-EMF  uB,  and  j  as  a 
function  of  distance  for  the  supersonic 
section  (approximately  latter  12  cm.)  of 
the  channel.  These  have  been  computed 
for  a  total  current  of  8030  A,  and  mass 
flow  of  2  g/s.  These  profiles  then 
provide  the  necessary  boundary 
conditions  for  solving  the  system  of 
equations  (3)  through  (10),  for  the 
cathode  temperature  distributions. 
Computations  using  the  method  outlined 
in  section  III  were  carried  out  on  a 
10x10  uniform  grid.  Finer  grids  (50x50 
and  100x100)  resulted  in  quantitatively 
different  temperatures,  but  trends  were 
identical.  Furthermore,  the  significant 
results  obtained  (see  Fig.  9)  were 
independent  of  the  grid  size.  As  can  be 
seen  from  Fig.  9,  the  maximum 
temperature  occurs  in  the  middle  of  the 
thruster  and  not  at  the  cathode  tip. 
Also,  the  internal  temperatures  are 
higher  than  the  surface  temperatures, 
due  to  significant  Joule  heating.  The 
surface,  once  molten  due  to  the  thermal 
runaway  caused  by  excessive  electron 
bombardment[7,8],  evaporates.  The 
interior  of  the  cathode,  on  the  other 
hand  once  molten  continues  to  change 
phase  and  build  up  pressure  within  the 
material.  This  explains  the  observed 
damage  to  the  ZT-1  and  DT-2  Stuttgart 
thrusters. 


V.  Back-EMF  Onset  and  Erosion 

Let  us  review  the  phenomena  occuring 
as  the  total  current  is  increased  for  a 
fixed  mass  flow  rate,  and  \  given 
geometry: 

(1)  The  back-EMF  climbs  from  a  small 
value  at  the  inlet,  reaches  a  maximum 
somewhere  in  the  middle  of  the  channel, 
and  then  decreases  toward  the  exit  (see 
Fig.  7). 

(2)  This  causes  the  net  current  density 
to  reach  a  minimum  in  the  middle  of  the 
channel,  while  being  large  near  exit 
region  (see  Fig.  8). 

(3)  Earlier  cathode  and  anode  sheath 
theories  [7,8]  have  shown  that  the 
diffuse  mode  behavior  becomes  unstable 
(possibly  transitioning  into  the  spot 
mode)  at  the  cathode  at  low  current 
densities  ,  and  at  the  anode  at  high 
current  densities. 

Thus,  the  center  of  the  cathode  and  the 
exit  regions  of  the  anode  are  prime 
candidates  for  regions  of  significant 
erosion. 

VL _ Back-EMF  Onset  and  Voltage 

Oscillations 

It  can  be  shown  from  current 
conservation,  that  the  sheath  voltage 
drop  follows  the  same  trend  as  the  net 
current  density[8].  The  regions  of  low 
current  density  at  the  cathode  then  have 
small  sheath  voltage  drops,  thereby 
making  the  thermal  runaway  more 
likely.  However,  the  exit  region  of  the 
cathode  has  higher  sheath  voltage  drops. 
It  is  therefore  clear,  that  as  the  total 
current  is  increased,  the  sheath  voltage 
drop  in  the  middle  of  the  cathode 
continuously  decreases,  while  the 
sheath  voltage  drops  near  the  exit 
region  increases.  When  the  sheath 
voltage  drop  reaches  a  value  comparable 
with  the  energy  of  the  first  excited 
state  of  the  propellant  atoms,  the 
electron  energy  distribution  function  is 
modified  near  the  electrode.  Two 
principal  groups  of  electrons  then 
result.  These  are  the  "low"  energy 
plasma  electrons,  and  the  energetic 
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electrons  emitted  by  the  electrode  and 
which  are  accelerated  through  the 
cathode  sheath.  Interaction  between 
these  two  groups  of  electrons  can  then 
lead  electric  field  or  voltage 

oscillations  arising  from  the  well 
known  beam  instability[18]. 

VII.  Summary 

A  method  has  been  outlined  that  allows 
determination  of  the  MPD  flow 
structure  and  electrode  temperature 

distributions.  given  the  global 

parameters,  mass  flow  rate.  Total 
current,  and  geometry.  This  theory 
extends  previous  work  on  back-EMF 
Onset,  to  include  erosion.  The  general 
predictions  of  this  theory  appear  to  be 
in  agreement  with  experimental 

onservations  at  Stuttgart[17], 
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Fig.  1:  Schematic  of  blown  cathode  from  DT-2 
thruster.  For  photograph,  see  ref.[17). 


Fig.  2:  Schematic  of  blown  cathode  from  ZT-1 
thruster  For  photograph,  see  ref. (17). 


Fig.  3:  The  calculated  values  of 
magnetic  induction  are  shown  here 
versus  distance  along  the  cathode,  for 
the  ZT-1  thruster. 


Fig.  4:  The  calculated  vafues  of 

velocity  are  shown  here  versus  distance 
along  the  cathode,  for  the  ZT-1 
thruster. 


Fig.  5:  The  calculated  values  of 

Temperature  are  shown  here  versus 

distance  along  the  cathode,  for  the  ZT-1 
thruster. 
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Fig.  6:  The  calculated  values  of 

ionization  fraction  are  shown  here 
versus  distance  along  the  cathode,  for 
the  ZT-1  thruster. 
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Fig.  7:  The  calculated  values  of  the 
back-EMF  (uB)  are  shown  here  versus 
distance  along  the  cathode,  for  the  2T-1 
thruster. 
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Fig.  8:  The  calculated  values  of  net 
current  density  are  shown  here  versus 
distance  along  the  cathode,  for  the  ZT-1 
thruster. 


Fig.  9:  The  calculated  centerline  and 
surface  temperatures  are  shown  here 
for  the  ZT-1  thruster.  For  reference, 
the  melting  temperature  of  tungsten 
(3680  K)  has  been  plotted  as  a 
horizontal  line.  Note  that  the 
temperatures  in  the  middle  region  are 
in  excess  of  the  melting  temperature, 
and  that  centerline  temperatures  are 
always  larger  than  the  surface 
temperature. 
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Electrode  erosion  in  magnetopiasmadynamic  thrusters  is  significantly  dependent  on  whether  the  mode  of 
operation  is  diffuse  (i.e.,  the  net  current  density  is  distributed  over  the  electrode  surface)  or  via  spots.  Since 
spot  mode  erosion  rates  are  typically  much  higher  than  the  diffuse  mode  erosion  rates,  the  latter  is  more  desirable 
from  the  standpoint  of  sustained  steady-state  operation.  In  this  paper,  the  steady-state  thermal  response  of  the 
cathode  in  a  magnetoplasmadynamic  thruster  is  analyzed  in  order  to  determine  the  limits  of  diffuse  mode 
operation.  An  energy  balance  at  the  cathode  surface  along  with  current  conservation  serve  to  simultaneously 
determine  the  surface  temperature  and  the  total  sheath  voltage  drop  at  a  given  location  on  the  electrode.  It  is 
found  that  two  different  limits  to  steady-state  diffuse  mode  operation  exist.  One  limit  corresponds  to  an  unsteady 
thermal  runaway  caused  by  excessive  electron  bombardment,  in  which  regenerative  heating  leads  to  local  melting 
of  the  electrode  surface.  The  other  limit  corresponds  to  an  operating  regime  where  a  steady  state  cannot  exist 
because  of  increased  ion  bombardment.  In  this  regime,  a  small  increase  in  the  net  plasma  current  density  results 
in  a  large  sheath  voltage  drop.  This  causes  an  increase  in  local  surface  heating  through  ion  bombardment,  until 
a  new  state  is  reached,  which  is  inherently  unsteady  and  accompanied  by  plasma  oscillations.  These  two  limits 
found  here  may  correspond  to  transitions  from  diffuse  to  prespot  modes  of  electrode  operation.  These  operating 
limits  are  also  found  to  be  influenced  strongly  by  electrode  and  discharge  parameters,  as  well  as  by  external 
cooling. 


Nomenclature 

A  =  thermionic  emission  constant 
C,  -  absolute  velocity  of  an  electron 
C„  =  component  of  the  absolute  velocity  of  an  electron  in 
the  y  direction,  i.e.,  in  the  direction  normal  to  the 
surface 

£  =  electric  field 

Ec  -  magnitude  of  the  v  component  of  the  electric  field 
at  the  cathode  surface 

£„  -  y  component  of  the  electric  field  at  the  sheath  edge 

e  =  electronic  charge 

f,  =  electron  velocity  distribution  function 

h  =  heat  transfer  coefficient  roughly  representing  overall 
heat  transfer  to  a’>  external  coolant  that  is  at  a  tem¬ 
perature  T coo, 

je  =  surface  electron  emission  current  density 
j,  =  current  density  of  plasma  electrons 

j,  =  current  density  of  plasma  ions 

/_  =  net  plasma  currer.  density  at  a  given  axial  location 

K,  =  thermal  conductivity  of  cathode  material 
k  =  Boltzmann's  constant 

L  =  cathode  length 

m,  =  mass  of  an  electron 

m,  =  mass  of  an  ion 

nt  =  electron  number  density 

n,  =  ion  number  density 
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nx  =  charged  particle  number  density  at  the  edge  of  the 
collisionless  sheath,  at  a  given  axial  location 
Q  =  net  heat  per  unit  area  per  unit  time  into  the  cathode 
surface,  at  a  given  axial  location 
T  =  local  cathode  surface  temperature 

Ta,„  =  temperature  of  an  external  coolant 
Tx  =  plasma  electron  temperature  at  a  given  axial  location 
V  =  voltage  or  potential 

V,  =  value  of  the  voltage  at  the  cathode  surface  at  a  given 
axial  location,  or  the  total  voltage  drop  across  the 
sheath 

y  =  coordinate  in  the  direction  normal  to  the  cathode 
surface 

f,  =  ionization  potential  of  the  propellant  gas 
<b  =  electrode  work  function 


I.  Introduction 

ELECTRODE  erosion  is  of  primary  importance  for  the 
prediction  of  lifetimes  of  magnetoplasmadynamic  (MPD) 
thrusters.  Erosion  processes  depend  on  a  complex  coupling 
between  plasma  discharge  characteristics,  plasma-wall  inter-’ 
actions,  and  electrode  phenomena.  In  particular,  erosion  rates 
depend  on  whether  the  current  conduction  is  through  spots, 
or  via  a  diffuse  (distributed)  mode.  The  diffuse  mode  is  char¬ 
acterized  by  distributed  current  emission  over  the  electrode 
surface,  surface  temperatures  well  below  the  material's  melt¬ 
ing  temperature,  and  non-negligible  plasma  ion  current.  By 
contrast,  spots  are  characterized  by  locally  constricted  or  fil¬ 
amentary  current  conduction,  surface  temperatures  at  or  above 
the  material's  melting  temperature,  and  strong  temperature- 
field  emission.  Spots  arc  detrimental  to  the  electrode  material 
because  of  their  high  erosion  rates.'  Therefore,  it  is  important 
to  understand  how  and  under  what  conditions  they  may  be 
formed  and  exactly  when  diffuse  mode  behavior  ends.  Current 
understanding  of  this  transition  from  diffuse  to  spot  behavior 
is.  at  best,  rather  poor.  In  this  paper,  the  steady-state  diffuse 
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mode  at  the  cathode  is  examined  theoretically.  Limits  on 
thruster  operation  in  this  mode  are  found,  and  an  explanation 
for  a  possible  transition  to  spot  formation  is  given.  Addition¬ 
ally,  the  analysis  yields  a  means  of  estimating  cathode  surface 
temperature,  which  can  be  used  to  predict  evaporative  erosion 
rates  for  the  cathode. 

Earlier  research  on  cathode  processes  has  focused  mainly 
on  spot  behavior  and  spot  characteristics.210  Although  these 
works  reveal  the  intricate  and  complex  phenomena  occurring 
in  a  cathode  spot,  little  or  no  information  is  provided  as  to 
how  spots  are  formed  in  the  first  place.  Some  authors  have 
attempted  to  describe  transition  from  diffuse  to  spot  modes. 
Moizhes  and  Rybakov"  have  found  a  negative-slope  region 
in  the  emission  current  voltage  characteristic,  which  they  at¬ 
tribute  to  transition  to  a  spot  mode.  These  authors  attribute 
the  transition  to  a  thermal  instability,  whose  physical  meaning 
or  origin  is  totally  unclear.  A  thermal  runaway  mechanism  in 
spots  due  to  a  positive  feedback  between  Joule  heating  and 
a  temperature  dependent  electrical  conductivity  has  been  pro¬ 
posed  by  Hantzsche.12  In  the  same  paper,  however,  the  author 
shows  that  thermionic  emission  cools  the  surface  and  prevents 
this  thermal  runaway.  It  is  important  to  note  that  Hantzsche 
assumes  the  sheath  voltage  drop  to  be  constant  and  given. 

By  contrast  to  the  cathode,  transition  at  the  anode  appears 
to  have  been  studied  more  extensively.11'11'  These  works  at¬ 
tribute  spot  formation  to  an  initial  local  surface  meltdown 
arising  primarily  from  Joule  heating  at  high  current  densities. 
However,  another  thermal  runaway  mechanism  involving  the 
much  ignored  sheath  exists.20  In  this  theory,  electron  emission 
from  the  electrode  surface  decreases  the  sheath  potential  drop, 
resulting  in  increased  electron  bombardment  from  the  plasma. 
This  results  in  a  positive  feedback  between  increased  emission 
and  heating  due  to  electron  bombardment,  causing  a  thermal 
runaway.  This  may  occur  when  the  anode  sheath  voltage  dif¬ 
ference  (defined  here  as  anode  potential  minus  the  adjacent 
plasma  potential)  is  negative  and  increasing.  This  thermal 
runaway  mechanism  arises  before  the  well-known  anode  sheath 
reversal21  and  may  explain  formation  of  anode  spots.  In  this 
paper,  the  same  mechanism  is  shown  to  be  also  operative  at 
the  cathode. 

A  simple  model  of  the  cathode  surface  and  its  adjacent 
sheath  is  discussed  in  Sec.  II.  Results  for  a  purely  thermion- 
ically  emitting  cathode  under  a  range  of  conditions  corre¬ 
sponding  to  MPD  thruster  operation,  are  given  in  Sec.  III. 
The  physical  meanings  of  the  limits  on  steady-state  diffuse 
mode  operation  are  discussed  in  Sec.  IV.  and  the  effects  of 
field-enhanced  and  temperature-field  emission  are  discussed 
in  Sec.  V.  Finally,  a  summary  of  this  work  along  with  its 
conclusions  are  presented  in  Sec.  VI. 


II.  Simple  Model 

This  section  will  focus  on  a  thermal  model  of  the  cathode 
surface  at  steady  state.  The  cathode  surface  is  subjected  to 
both  heating  and  cooling  (Fig.  1).  Charged  particles  (i.e.. 
electrons  and  ions)  from  the  plasma  bombard  the  surface  and 
consequently  heat  it.  The  ions  (which  we  will  assume  to  be 
singly  charged)  also  recombine  with  electrons  at  the  electrode . 
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Fig.  I  Cathode  along  with  the  various  heating  and  cooling  mecha¬ 
nisms. 


thereby  heating  the  surface  even  further.  Additional  heating 
due  to  radiation  from  the  plasma  is  also  present  and  expected 
to  be  important,22  but  will  be  neglected  here.  The  electrode 
may  also  emit  electrons,  which  can  result  in  heating  (during 
field  emission)  or  in  cooling  (during  thermionic  emission). 
During  steady-state  operation,  it  is  likely  that  the  electrode 
is  cooled  externally.1 23  This  is  usually  done  by  cooling  the 
cathode  at  its  base.  Finally,  the  cathode  surface  may  radiate 
away  its  heat.  With  these  considerations,  we  may  proceed  to 
write  the  following  for  the  net  heat  area  per  unit  and  per  unit 
time  entering  the  cathode  surface: 

Q  =  i.  [«*>  +  +  J.(K  +  «.-*) 

-  i,  [d>  +  v"]  -  *<r  ~  r-,)  (i) 

The  quantity  h(T  -  Tcool)  represents  a  very  rough  model  of 
the  overall  cooling,  and  h  can  be  considered  as  equivalent  to 
the  electrode  material  thermal  conductivity  per  unit  length 
(i.e.,  KJL).  At  steady  state,  Q  must  be  0  and  Eq.  (1)  deter¬ 
mines  the  surface  temperature  T.  Radiative  heat  transfer  to 
the  anode  can  be  shown  to  be  negligible  and  is.  therefore, 
ignored  here. 

Consider  next  the  sheath  region  immediately  adjacent  to 
the  electrode  surface.  This  region  is  typically  of  the  order  of 
tens  of  Debye  lengths.  The  electron  mean  free  path,  however, 
is  orders  of  magnitude  larger.22  Consequently,  the  sheath  may 
be  treated  as  a  collisionless  sheath.  The  electrons  entering  the 
sheath  from  the  plasma  are  then  described  by  the  collisionless 
steady  Boltzmann  equation,  which  in  one  dimension  is 


Using  the  following  transformations. 


,  mtC\  ,  eV 
-  2kT,  kT 


t\  =  y 


(2) 


(3) 


along  with  E  =  —HViby,  yields 


—  =  0.  or  that  /,  =  /,(£)  onlv  (4) 


With  the  condition  that  the  electron  distribution  function  be 
Maxwellian  at  the  sheath  edge  (where  V  is  taken  to  be  0).  we 
obtain  the  following  exact  solution  for  the  electron  distribu¬ 
tion  function  in  the  sheath: 


/,=  nr(mr!2nkTrY  2  exp(  —  m,C2'2fc7\].  (5) 

where 

n,  =  n,  exp[  -  eV/kT,]  (6) 

with  n ,  being  the  electron  number  density  at  the  edge  of  the 
sheath  where  V  is  taken  to  be  0.  The  average  flux  of  electrons 
may  then  be  calculated  easily  from  Eqs  (5)  and  (6)  by  in¬ 
tegrating  fr  over  dC,.y  and  dCr.  from  -  *  to  *.  and  over  dCry 
from  0  to  *.  This  flux  when  evaluated  at  the  wall  and  mul¬ 
tiplied  by  e  gives  the  electron  current  density  at  the  wall: 

/ ,  =  en,{kTj2m,ir)' 2  exp[  -  eV,  kTx]  (7) 

where  -  V,  is  the  potential  of  the  cathode  with  respect  to  the 
sheath  edge  (where  V  =  0).  and  /,  is  taken  to  be  positive  in 
the  direction  away  from  the  cathode. 

The  ions  entering  the  sheath  are  assumed  to  be  mono- 
energetic.  The  Bohm  criterion  for  a  stable  monotonic  sheath 
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is,  therefore,  satisfied,24  This  condition  is  only  slightly  altered 
in  the  presence  of  electron  emission  from  the  surface.25  More 
exact  theories  that  provide  the  velocity  distribution  function 
of  the  ions  entering  the  sheath  from  the  cotlisional  presheath 
exist.26-2*  However,  for  present  purposes,  the  frequently  em¬ 
ployed  monoenergetic  assumption  simplifies  the  problem  a 
great  deal.  We  may  then  write  the  ion  current  density  at  the 
surface  as 


j,  =  enJLkTJm.Y*  (8) 

where  j,  is  taken  to  be  positive  toward  the  cathode  and  quasi¬ 
neutrality  ( n ,  =  nr  =  n.)  is  assumed  as  the  sheath  edge.  The 
sheath  edge  is  a  rather  arbitrary  definition  since  no  such  rigid 
demarcation  exists  in  reality.  When  the  ions  are  monoener¬ 
getic,  the  sheath  edge  has  a  clear  meaning.  However,  even 
when  the  ion  distribution  function  is  not  monoenergetic, 
the  sheath  edge  can  be  defined  as  the  location  where  the 
velocity  of  the  ions  reaches  the  value  of  the  Bohm  velocity 
(kr*/m,)1/2.  In  this  paper,  we  will  adopt  this  definition  and 
take  the  data  for  the  sheath  potential  at  this  location  (i.e.,  V 
=  0  at  the  sheath  edge).  Now,  if  the  surface  emission  current 
density  denoted  by  jc  is  taken  to  be  positive  toward  the  surface 
and  the  net  plasma  current  density  denoted  by  /„  is  taken  to 
be  positive  toward  the  cathode,  we  may  write  overall  current 
conservation  at  steady  state  as 


A  =  )e  +  i.  ~  A  (9) 

where,  in  general,  jE  is  a  function  of  the  local  surface  tem¬ 
perature  T  and  the  local  electric  field  at  the  cathode  surface 
Ec.  Determination  of  Ec  is  discussed  in  Sec.  V.  Equation  (9) 
determines  the  total  cathode  sheath  voltage  drop  Vc.  Given 
the  parameters  7"x.  </>,  h,  Tcool  and  the  surface  electric 

field  Ec.  Eqs.  (1)  and  (9)  represent  two  simultaneous  equa¬ 
tions  in  the  two  unknowns  T  and  Ve.  It  must  be  pointed  out 
that  Crawford  and  Cannara  mention  overall  current  conser¬ 
vation  in  their  work.2v  However,  their  regime  of  interest  was 
at  very  low  densities  ( 1014  m  '),  and,  hence,  ion  and  electron 
currents  from  the  plasma  were  negligible  in  comparison  with 
the  emission  current.  This  is  an  important  difference  between 
their  work  and  the  present  work. 

A  simple  solution  may  be  found  for  the  case  of  pure  therm¬ 
ionic  emission  (which  is  likely  for  refractory  materials),  where 
the  surface  emission  current  density  depends  only  on  T.  For 
this  case,  the  total  sheath  voltage  drop  may  be  explicitly  ob¬ 
tained  from  Eq.  (9): 


(10) 


where  /,  =  en%(kTJ2itm,)' :.  Substituting  Eq.  (10)  into  Eq. 
(1)  yields  the  following: 


Q  = 


T)  +  /, 


+  e,  + 


-  A 


-  h(T  -  7_t) 


(11) 


Since  at  steady  state  Q  must  be  zero,  Eq.  (11)  represents  an 
implicit  equation  for  the  local  surface  temperature  T,  where 
V,  is  given  by  Eq.  (10). 

An  interesting  phenomenon  is  revealed  by  Eq.  (11).  The 
emission  current  density,  which  under  thermionic  conditions 
results  in  energy  transport  away  from  the  cathode,  appears 
as  an  input  heating  source.  This  can  be  understood  upon  close 
examination  of  Eq.  (9)  or  (10).  When  the  emission  current 
density  increases  for  a  fixed  plasma  number  density  and  plasma 
current  density,  the  sheath  voltage  drop  must  decrease  in 
order  to  conserve  total  current.  This  results  in  an  increase  in 


electron  bombardment,  which  contributes  to  heating  the  sur¬ 
face.  Thus,  although  thermionic  emission  cools  the  surface, 
its  overall  effect  is  to  heat  the  surface  via  increased  electron 
bombardment.  It  is,  therefore,  important  to  include  the  ef¬ 
fects  of  the  sheath  when  attempting  to  determine  the  surface 
temperature. 

A  simple  model  has  been  developed  in  this  section,  which 
considers  a  thermal  balance  of  the  cathode  surface  along  with 
overall  current  conservation.  The  local  surface  temperature 
and  overall  sheath  voltage  drop  can  be  determined  from  this. 
Results  for  a  purely  thermionically  emitting  cathode  under 
typical  MPD  conditions  are  presented  next. 

III.  Results  for  Pure  Thermionic  Emission 

The  governing  equations  presented  in  the  preceding  section 
are  solved  in  this  section  for  the  case  of  pure  thermionic 
emission  (i.e.,  no  field-enhanced  or  temperature-field  emis¬ 
sion).  Temperatures  and  sheath  voltage  drops  are  calculated 
for  typical  MPD  conditions. 

For  pure  thermionic  emission,  the  emission  current  density 
is  given  by 


jE  =  AT2  exp{-e<t>/kT]  (12) 

The  sheath  voltage  drop  Ve  can  then  be  written  in  terms  of 
T  using  Eq.  (10).  Given  the  parameters  wx.  yx,  7"x,  <J>,  h,  and 
rcool,  Eqs.  ( 10—12)  may  be  combined  and  solved  for  T.  The 
sheath  voltage  drop  can  then  be  determined  from  Eq.  (10). 
The  emission  coefficient  may  vary  considerably  for  materials 
with  oxide  layers.  Both  A  and  <j>  may  vary  locally  on  a  given 
surface,  so  that  it  is  important  to  consider  the  sensitivity  of 
final  results  to  variations  in  these  parameters.  This  is  ad¬ 
dressed  later  in  this  section. 

A  typical  variation  of  the  net  heat  into  the  cathode  Q  vs 
the  surface  temperature  is  displayed  in  Fig.  2.  Two  intersec¬ 
tions  with  the  horizontal  axis  are  found,  representing  two 
possible  steady-state  solutions.  The  first  (lower  temperature) 
is  a  stable  attractor  and  the  second  (higher  temperature)  is 
an  unstable  repeller.  This  can  be  seen  quite  easily  by  per¬ 
turbing  the  solutions  to  either  side  and  determining  if  the 
initial  state  is  restored.  The  physical  meaning  of  the  stable 
point  is  clear.  The  incoming  energy  on  the  surface  is  exactly 
balanced  by  the  outgoing  energy.  Furthermore,  the  surface 
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Fig.  1  Typical  variation  of  the  net  heal  per  unit  area  per  unit  time 
into  the  cathode  surface  vs  surface  temperature.  Variation  predicted 
by  Eqs.  (10)  and  (II)  in  Sec.  II. 
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can  maintain  this  temperature  for  an  indefinite  period  of  time, 
even  for  small  fluctuations  in  the  discharge.  The  unstable 
point  requires  some  explanation.  Although  an  exact  balance 
between  incoming  and  outgoing  energy  is  possible  at  this  point, 
the  surface  cannot  remain  in  this  state  indefinitely  and  is 
susceptible  to  change  due  to  infinitesimal  fluctuations.  Phys¬ 
ically,  this  state  represents  an  operating  regime  where  the 
temperature  and,  hence,  the  emission  current  density  is  high 
enough  to  lower  the  sheath  voltage  drop  below  a  critical  value. 
This  causes  excessive  heating  due  to  increased  electron  bom¬ 
bardment,  resulting  in  a  temperature  rise.  The  subsequent 
increase  in  temperature  causes  an  increase  in  current  emis¬ 
sion,  which  lowers  the  sheath  voltage  drop  further.  This  pos¬ 
itive  feedback  process  then  repeats  itself  until  the  surface  is 
regeneratively  heated  up  to  the  melting  temperature.  This 
effect  is  quantitatively  described  in  Sec.  IV. 

Results  from  the  present  calculation  are  shown  in  Figs.  3- 
6.  In  these  figures,  the  steady-state  surface  temperatures  cal¬ 
culated  from  Eq.  (11)  are  shown  vs  the  net  plasma  current 
density  /«,  the  charged  particle  number  density  at  the  sheath 
edge  n_,  the  electrode  work  function  <f>,  and  the  temperature 
of  the  external  coolant  The  heat  transfer  coefficient  h 
is  fixed  at  20  kW/mJ/K,  the  plasma  electron  temperature  is 
taken  to  be  12,000  K,  and  A  =  3  x  10*  A/m2/K2.  The  value 
of  h  is  obtained  from  an  estimate  of  the  thermal  conductance 
through  a  tungsten  cathode.  All  four  plots  display  both  the 
stable  (lower  temperature)  and  unstable  (higher  temperature) 
regimes.  Figures  3  and  4  display  the  varition  of  T  with  and 
ft.,  respectively,  for  <t>  =  2.63  V  and  =  500  K.  From 
Fig.  3,  it  can  be  seen  that,  for  a  given  number  density,  the 
stable  diffuse  mode  is  limited  both  at  low  as  well  as  high  values 
of  the  current  density.  Similarly,  for  a  given  current  density, 
there  are  lower  and  upper  limits  on  the  range  of  allowable 
number  densities.  At  low  current  densities  (n„  fixed)  and  at 
high  number  densities  (/„  fixed),  the  stable  and  unstable  so¬ 
lutions  are  seen  to  merge  and  disappear.  Also,  at  the  higher 
current  densities  ( « .  fixed)  and  the  lower  number  densities 
(/'„  fixed),  the  stable  diffuse  operation  is  limited  and  only  the 
unstable  mode  exists.  It  is  interesting  to  note  that,  for  a  given 
number  density,  the  steady-state,  stable  surface  temperature 
does  not  vary  significantly  with  the  current  density  over  a 
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Fig.  3  Steady-state  surface  temperature  as  predicted  by  Eqs.  (10) 
and  (II)  vs  net  plasma  current  density  y,  for  various  values  of  n,:  h 
=  20  kW/mJ/K;  *  =  2.63  V;  Tc„<  =  500  K;  A  =  3  x  10*  A/m’/K'; 
r.  *  12000  K. 
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Fig.  4  Steady-state  surface  temperature  as  predicted  by  Eqs.  (10) 
and  (11)  vs  net  plasma  current  density  n»  for  various  values  of  j,s  h 
=  20  kW/m2/K;  *  =  2.63  V;  TM  *  500  KM  =  3  x  10*  A/mVK2; 
r.  =  12000  K. 


o 
o 
1 0 
n 


o 

u  n 

a 

4J 

2  o 
a  cv 


H  § 
o 
v  N 
0 

43  . 

1  ° 

h  ° 

3  «o 
cn  - 

v  . 
v  o 

2  o  h 


43 

o 


A=lx\(fAlmziK2 


// 

>-4=I20xlO,A/m2/K2 


Oxio11 


6  v 


qxio”w-> 


STABLE 


OO  to  8  0  3  0  4.0 

Cathode  Material  Work  Function 


5  O 


[V] 


Fig.  5  Steady-state  surface  temperature  as  predicted  by  Eqs.'  (10) 
and  (11)  vs  electrode  work  function  for  jm  =  10*  A/m1,  n„  =  8  x 
1021  m'J,  =  500  K,  A  =  20  kW/mVK,  and  F.  =  12000  K. 


wide  range  (see  Fig.  3).  This  is  because  in  this  region  the  ion 
current  primarily  determines  the  sheath  voltage  drop,  which 
is  then  to  a  very  good  approximation  a  constant  since  the 
number  density  is  fixed.  Heating  due  to  ion  bombardment  is 
then  balanced  by  the  external  cooling.  Therefore,  inspection 
of  Eq.  (11)  reveals  that  the  surface  temperature  should  not 
depend  significantly  on  /.  in  this  region.  By  contrast,  changes 
in  ri.  at  a  fixed  /,  affect  the  stable  operation  drastically  via 
increased  electron  and/or  ion  bombardment.  This  is  reflected 
in  Fig.  4.  The  variation  of  the  surface  temperatures  with  elec¬ 
trode  work  function  is  shown  in  Fig.  5  for  Tcool  =  500  K,  /. 
=  106  A/m2,  and  «.  =  8  x  1051  m~\  The  previously  discussed 
stable  and  unstable  regimes  are  again  found  to  exist.  The 
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plasma  current  density  and  plasma  number  density  at  the 
sheath  edge,  Eq.  (10)  may  be  substituted  into  Eq.  (11)  and 
the  resulting  expression  differentiated  with  respect  to  7  to 
yield 


(13) 


For  stability  (in  the  Lyapunov  sense),  we  require  that  dQ/dT 
be  less  than  zero.  A  sufficient  condition  for  stability  from  Eq. 
(13)  is 


(14) 


which,  by  using  Eqs.  (7)  and  (8),  gives  a  condition  on  the 
sheath  voltage  drop  Vc: 


Fig.  6  Steady-state  surface  temperature  as  predicted  by  Eqs.  (10) 
and  (11)  vs  external  coolant  temperature  for  <6  =  2.63  V,J.  = 
10*  A/mJ,  h  =  20  kW/nrVK,  A  =  3  x  10*  A/mVKJ,  and  T.  =  12000 
K. 


stable,  steady-state  surface  temperature  is  almost  invariant 
with  respect  to  <f>  in  this  regime,  again  because  the  ion  current 
principally  determines  the  sheath  drop  and  ion  bombardment 
dominates.  Also  shown  in  this  figure  is  the  solution  for  a 
different  value  of  the  thermionic  emission  constant  (A  =  120 
x  10*  A/m2/K2).  The  effect  of  a  higher  value  of  A  is  simply 
to  shift  the  turning  point  (i.e.,  the  point  where  stable  and 
unstable  solutions  merge  and  then  disappear)  toward  the  higher 
work  function  side.  Otherwise,  the  stable,  steady-state  surface 
temperature  is  practically  the  same  as  for  a  lower  value  of  A 
since  ion  bombardment  is  still  the  dominant  mechanism.  Fig¬ 
ure  6  displays  the  surface  temperatures  vs  the  external  coolant 
temperature  for  <t>  =  2.63  V,  /„  =  106  A/m2,  nm  -  5  x  1021 
m'3,  and  h  =  20  kW/m2/K.  Both  stable  and  unstable  regimes 
are  seen  to  exist  and  exhibit  the  features  that  have  just  been 
discussed. 

This  section  has  focused  on  results  for  a  thermionically 
emitting  cathode  operating  under  typical  MPD  conditions. 
Stable  and  unstable  regimes  as  well  as  limits  on  steady-state 
operation  have  been  found.  These  limits  are  found  to  be 
influenced  strongly  by  both  discharge  and  electrode  param¬ 
eters.  The  next  section  will  focus  on  a  detailed  discussion  of 
the  operating  limits  that  have  been  found. 


IV.  Limits  on  the  Steady  Diffuse  Mode 
In  this  section,  two  different  mechanisms  that  lead  to  local 
melting  of  the  cathode  surface  are  discussed.  The  first  mech¬ 
anism  is  an  excessive  heating  due  to  electron  bombardment 
that  occurs  when  the  sheath  voltage  drop  falls  below  a  critical 
value.  A  positive  feedback  between  electron  bombardment 
and  surface  electron  emission  results  in  this  thermal  runaway. 
The  second  mechanism  is  heating  due  to  ion  bombardment, 
which  dominates  as  the  sheath  voltage  drop  increases  sharply 
for  a  relatively  small  increase  in  the  plasma  current  density. 
A  steady  state  can  no  longer  be  maintained  unless  plasma 
discharge  parameters  change  dramatically. 

We  consider  the  thermal  runaway  mechanism  first.  The 
value  of  the  critical  voltage  drop  at  which  the  thermal  runaway 
due  to  electron  bombardment  occurs  can  be  obtained  quite 
readily  for  the  case  of  pure  thermionic  emission.  For  a  given 


For  argon  and  an  electron  temperature  of  12,000  K,  this  yields 
an  approximate  value  of  5.5  V  for  the  critical  sheath  voltage 
drop.  It  is  important  to  note  that  this  is  only  a  sufficient 
condition  for  stability.  From  Eq.  (13),  it  can  be  seen  that  the 
exact  value  of  the  critical  voltage  drop  depends  on  the  elec¬ 
trode  properties  (work  function  and  emission  current),  sur¬ 
face  temperature,  and  external  cooling  h.  This  more  general 
criterion  on  the  total  sheath  voltage  drop  necessary  to  ensure 
stability  of  a  given  steady  state  is 


7(3  -t-  e<t>/kT) 
T„(2  +  e<t>/kT) 


_ ehT _ \ 

kjET^(2  +  e<j>/kT)/ 


(16) 


Above  the  value  given  by  the  critical  voltage  drop,  a  steady 
state  defined  by  Eqs.  (10)  and  (11)  (with  Q  =  0)  can  be 
maintained  indefinitely,  even  for  small  fluctuations  in  and 
/«.  Below  this  value,  the  positive  feedback  between  electron 
bombardment  and  surface  electron  emission  leads  to  thermal 
runaway,  and  a  stable  steady  state  cannot  be  maintained.  The 
surface  subsequently  melts  locally.  Alternatively,  the  condi¬ 
tion  dQIBT  <  0  may  be  interpreted  as  providing  a  constraint 
on  the  amount  of  cooling  (i.e.,  the  value  of  h)  necessary  for 
steady  and  stable  operation.  Equation  (14)  along  with  Eq.  (9) 
may  be  rewritten,  therefore,  in  the  following  form,  yielding 
an  important  upper  limit  on  the  local  surface  temperature  for 
diffuse  operation: 

ic(T)  =  AT-  exp(-  e<t>lkT)  s  jm  -  j,  (17)  ' 

Under  MPD  conditions,  the  critical  temperatures  predicted 
by  Eq.  (17)  are  well  below  the  material  melting  temperature. 

A  second  mechanism  exists  that  can  cause  local  melting  of 
the  cathode  surface.  For  high /„  (given  an  nj)  or  low  n,  (given 
a  /'.),  it  is  possible  that  -*  j,  +  jE.  When  this  occurs,  the 
total  sheath  voltage  drop  increases  sharply.  This  can  be  seen 
from  Eq.  (10).  where  Vc  -»  *  as  /«  — *  j,  +  jE.  This  has  two 
important  consequences.  First,  emitted  electrons  cause  in¬ 
creased  ionization  outside  the  sheath  since  they  gain  large 
amounts  of  energy  after  traversing  the  sheath.  Second,  the 
ions  gain  substantial  amounts  of  energy  while  falling  through 
the  sheath.  These  cause  a  rapid  increase  in  ion  bombardment 
with  subsequent  heating  of  the  surface.  This  temperature  rise 
necessarily  results  in  increased  surface  electron  emission. 
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Consequently,  there  arises  a  conflict  between  the  emission 
current  necessary  for  cooling  the  cathode  surface  (i.e.,  to 
balance  heating  by  ion  bombardment)  and  the  emission  cur¬ 
rent  necessary  to  preserve  overall  current  conservation  at  steady 
state.  The  increasing  sheath  drop  also  drastically  reduces  the 
electron  current  from  the  plasma  necessary  to  counteract  an 
increasing  emission  current  and,  hence,  maintain  a  steady 
state.  This  leads  to  unsteady  behavior  in  the  sheath,  thereby 
forcing  the  surface  and  the  plasma  to  also  become  unsteady. 
In  the  meanwhile,  the  local  surface  temperature  can  rise  due 
to  ion  bombardment  until  the  surface  has  locally  melted.  Al¬ 
though  the  present  theory  becomes  inapplicable  at  this  point, 
it  is  capable  of  predicting  this  limit.  An  important  phenom¬ 
enon  associated  with  this  second  mechanism,  but  not  with  the 
first,  is  the  triggering  of  plasma  oscillations.30  When  the  sheath 
voltage  drop  rises  and  approaches  energies  equal  to  or  greater 
than  that  of  the  first  excited  states  of  the  propellant  atoms, 
ionization  levels  can  be  enhanced  in  the  presheath  region. 
Also,  while  the  thermionically  emitted  electrons  relax  their 
momenta  relatively  quickly  due  to  elastic  collisions  with  neu¬ 
trals,  their  energy  relaxation  times  are  much  longer  due  to 
infrequent  collisions  with  plasma  electrons.  When  the  number 
density  of  thermionically  emitted  (beam)  electrons  ap¬ 
proaches  the  local  number  density  of  plasma  electrons,  lon¬ 
gitudinal  oscillations  can  result.30  This  condition  may  be  ex¬ 
pressed  quantitatively  within  the  context  of  the  present  work 
as 

n,  «  2.563  x  1013  /g*Vg*  (18) 


o 


Heat  Transfer  Coefficient  H  (W/m*/fc) 


Fig.  8  Variation  of  the  cathode  surface  temperature  as  predicted  by 
Eqs.  (10)  and  (11)  vs  the  heat  transfer  coefficient  h  for  a  net  current 
density  ofj»  =  10*  A/m2  and  for  two  different  number  densities:  r 
=  12,000  K;  *  =  2.63  V;  =  500  K. 


where  n ,  is  the  number  density  of  plasma  electrons  immedi¬ 
ately  adjacent  to  the  sheath.  This  second  mechanism,  which 
is  related  to  a  high  back-EMF  (Electro  Motor  Force),  may 
explain  the  onset  phenomenon  observed  in  the  MPD  thrus¬ 
ters.3133 

The  two  aforementioned  mechanisms  responsible  for  lim¬ 
iting  the  diffuse  mode  also  cause  local  melting  of  the  cathode 
surface.  Local  surface  melting  may  be  a  precursor  to  spot 


formation.  However,  it  must  be  pointed  out  that  an  inter¬ 
mediate  stable  state  may  exist  where  molten  regions  that  are 
submicron  in  size  (microspots1)  are  scattered  over  the  cathode 
surface  while  most  of  the  discharge  appears  diffuse .  It  is  there¬ 
fore  possible  that  the  two  mechanisms  found  and  discussed 
here  correspond  to  transitions  from  diffuse  to  spot  mode  or 
another  intermediate  mode  with  microspots. 
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Fig.  7  Variation  of  the  electric  field  in  the  steady,  collisionless  elec¬ 
trode-adjacent  sheath  vs  vertical  distance  from  the  cathode  surface: 
/.  *  10*  A/m2;  n,  *  S  x  102'  m  J;  T,  =  12000  K;  Tc<^  =  500  K; 
♦  *  2.63  V;  h  =  20  kW/m2/K2;  the  calculated  surface  temperature 
is  1824  K;  the  value  of  the  electric  field  at  the  cathode  surface  is  1.76 
x  to2  V/m.  (Note  that  the  electric  field  is  negative  because  of  the 
convention  used  here.  The  actual  field  is  positive,  i.e.,  pointing  toward 
the  surface.) 


V.  Effects  of  Field-Enhanced  and 
Temperature-Field  Emission 

Thus  far,  we  have  considered  a  purely  thermionically  emit¬ 
ting  cathode.  In  general,  however,  the  emission  current  den¬ 
sity  depends  both  on  the  surface  temperature  and  on  the 
surface  electric  field.33  In  this  section,  the  influence  of  this 
field-enhanced  (also  known  as  Schottky  emission)  and  tem¬ 
perature-field  (T-F)  emission  is  evaluated  and  discussed. 

The  governing  equations  for  the  sheath  and  the  cathode 
surface  are  described  by  the  surface  energy  balance  [Eq.  (1)] 
and  overall  current  conservation  [Eq.  (9)j.  These  equations 
are  still  applicable  in  the  presence  of  an  electric  field  at  the 
surface  of  the  cathode.  However,  a  simple  analytic  solution 
cannot  be  obtained  because  the  sheath  voltage  drop  cannot 
be  determined  explicitly  in  terms  of  the  surface  temperature. 
This  is  because  the  emission  current  density  in  Eq.  (9)  depends 
on  T  and  Ec,  and  Er,  in  turn,  depends  on  Vc,  The  required 
additional  equation  for  Ec  may  be  obtained  from  a  solution 
of  Poisson’s  equation  in  the  sheath.  An  outline  of  this  pro¬ 
cedure  is  given  by  Prewett  and  Allen33  for  the  case  of  a  con¬ 
stant  surface  emission  current  density.  Extension  of  their  so¬ 
lution  for  the  potential  distribution  in  the  sheath,  for  the  case 
where  the  emission  current  density  depends  on  both  T  and 
Er.  is  straightforward  and  will  not  be  discussed  here.  The 
additional  equation  for  Ec  for  monoenergetic  ions  is  the  fol¬ 
lowing  implicit  equation: 


-  2  +  exp[-eV./*7\)j  -  jr(2m,  V,/ey  '- 


(19) 
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where  Ew  is  the  electric  field  at  the  sheath  edge  where  V  - 
0.  In  general,  Ew  is  much  smaller  than  Ee  and  may  be  ignored, 
and  the  emission  current  density  jE  may  be  calculated  from 
quantum-mechanical  considerations.  A  general  integral 
expression  for  jE  as  a  function  of  T  and  Ec  is  given  by  Murphy 
and  Good.'1  For  low  values  of  Ec,j£  may  be  simply  expressed 
by  the  well-known  Schottky  formula.3’  For  the  higher  values 
of  Ec,  the  emission  is  influenced  by  both  Ec  and  T  and  is 
known  as  T-F  emission.  Thus,  Eqs.  (1).  with  (2  =  0.  (9),  and 
(19)  represent  three  coupled  and  nonlinear  equations  for  the 
three  unknowns  T.  Vc,  and  Ec.  This  system  has  been  solved 
numerically.  It  is  also  possible  to  numerically  integrate  Pois¬ 
son’s  equation  to  determine  the  variation  of  the  potential  and 
the  electric  field  in  the  sheath.  Figure  7  displays  the  computed 
variation  of  the  electric  field  vs  vertical  distance  from  the 
surface  for  a  particular  value  of  ;‘x  and  n*.  However,  it  must 
be  pointed  out  that  the  resulting  profiles  represent  the  solution 
within  the  sheath  only.  A  uniformly  valid  solution  for  the 
sheath,  transition  region,  and  the  plasma  beyond  has  to  be 
obtained  through  matched  asymptotic  expansions. Con¬ 
sequently,  the  solution  of  Poisson’s  equation  given  here  rep¬ 
resents  only  the  inner  solution.  For  the  range  of  parameters 
considered  in  this  paper,  computed  values  of  E,  are  of  the 
order  of  107  V/m  or  less.  Such  values  of  the  electric  field 
produce  a  negligible  increase  in  current  emission  at  surface 
temperatures  of  interest  to  the  MPD  thruster.  Consequently, 
the  stable  diffuse  mode  is  well  represented  by  pure  thermionic 
emission. 

VI.  Discussion  and  Conclusions 

The  thermal  response  of  the  cathode  in  an  MPD  thruster 
has  been  studied  under  steady-state  diffuse  conditions.  The 
electrode-adjacent  sheath  has  been  included  for  monoener- 
getic  ions.  The  solution  of  Poisson’s  equation  in  the  sheath, 
together  with  overall  current  conservation  and  an  energy  bal¬ 
ance  at  the  electrode  surface,  serve  to  determine  the  electric 
field  at  the  surface,  the  surface  temperature,  and  the  total 
sheath  voltage  drop  simultaneously  and  self-consistently.  Two 
operating  regimes  are  found.  One  is  stable,  while  the  other 
is  an  unstable  thermal  runaway  caused  by  excessive  electron 
bombardment.  The  thermal  runaway  leads  to  eventual  local 
melting  of  the  cathode  surface  and  may  be  a  precursor  to  spot 
formation.  The  stable  steady  state,  on  the  other  hand,  is  in¬ 
fluenced  strongly  by  discharge  and  electrode  parameters.  The 
values  of  and  T,  used  here  have  been  obtained  from 

an  approximate  nonequilibrium  model  of  the  electrode-ad¬ 
jacent  flow.11  The  rate  of  cooling  is  found  also  to  be  extremely 
important.  Varying  values  of  h  result  in  two  operating  re¬ 
gimes,  similar  to  the  results  displayed  in  Fig  6.  Furthermore, 
the  dependence  of  the  surface  temperature  on  h  is  quite  sig¬ 
nificant.  For  =  5  x  1011  m'3,/x  =  10®  A/m:.  Tcool  =  500 
K,  and  <f>  =  2.63  V,  a  l %  decrease  in  h  results  in  about  a 
10%  increase  in  the  local  surface  temperature  (see  Fig.  8). 
Of  course,  for  h  larger  (=  105  W/m2/K)  than  the  value  used 
here,  the  variation  is  less  dramatic  but  still  non-negligible. 

Two  mechanisms  are  found  to  be  responsible  for  destabi¬ 
lizing  the  stable,  steady,  diffuse  mode.  The  first  occurs  at 
values  of  the  sheath  voltage  drop  below  a  critical  value.  Small 
sheath  voltage  drops  lead  to  increased  electron  bombardment 
and  subsequent  thermal  runaway  if  the  cooling  is  insufficient. 
Such  a  condition  can  occur  in  the  MPD  thruster  operating  at 
near-onset  total  currents.  Near  onset,  the  current  density  is 
sharply  peaked  and  large  near  the  inlet  and  the  exit  regions 
while  going  through  a  minimum  in  the  middle  region  (in  the 
axial  direction)  of  the  thruster.  The  current  density  in  the 
middle  region  of  the  thruster  becomes  smaller  with  increasing 
total  current  to  the  thruster,  due  to  a  high  back-EMF.31  '-’ 
which  then  leads  to  smaller  sheath  voltage  drops.  The  second 
mechanism,  due  to  excessive  ion  bombardment  and  accom¬ 
panied  by  plasma  oscillations,  occurs  at  high  values  of  the 
current  density  (for  a  given  plasma  number  density)  or  for 


low  values  of  the  plasma  number  density  (for  a  given  current 
density).  Here,  high  sheath  voltage  drops  result  in  increased 
ion  bombardment  leading  to  increased  surface  electron  emis¬ 
sion.  Since  plasma  electrons  are  repelled  at  high  sheath  po¬ 
tentials.  overall  current  conservation  given  by  Eq.  (9)  cannot 
be  satisfied  at  steady  state.  The  sheath  is  then  no  longer 
steady.  The  second  mechanism  may  be  inhibited  if  an  increase 
in  the  sheath  edge  charged  particle  number  density  n„  occurs 
via  increased  ionization  due  to  electrons  emitted  from  the 
cathode  surface.  Clearly,  this  requires  overall  sheath  voltage 
drops  of  the  order  of  the  propellant  ionization  potential.  The 
occurrence  of  both  mechanisms  is  strongly  influenced  by  dis¬ 
charge  parameters. 

The  simple  model  discussed  in  this  paper  represents  a  first 
attempt  at  relating  the  plasma  discharge,  the  sheath,  and  the 
electrode.  Although  this  model  has  revealed  some  important 
underlying  physics,  there  are  several  deficiencies.  First,  the 
presheath  has  not  been  considered  in  detail.  The  presheath 
is  important  not  only  for  determining  the  ion  velocity  distri¬ 
bution  function,  but  also  for  asymptotically  determining  a 
uniformly  valid  potential  distribution.  This  potential  distri¬ 
bution  will  then  yield  the  correct  variation  in  the  sheath  when 
considered  on  the  scale  of  the  Debye  length  and  will  also  yield 
the  correct  variation  in  the  plasma  when  viewed  on  the  mac¬ 
roscopic  length  scale  in  the  problem.  Such  presheath  solutions 
exist,  as  mentioned  previously  in  Sec.  2,  but  only  for  charge- 
exchange  reactions  or  in  the  absence  of  surface  electron  emis¬ 
sion.  For  the  plasma  in  the  MPD  thruster,  additional  ioni¬ 
zation  by  electrons  emitted  from  the  surface  is  important  under 
some  conditions.  Second,  radiative  heating  of  the  electrode 
surface  by  the  plasma  has  been  neglected.  This  not  only  con¬ 
tributes  to  direct  heating,  but  also  to  increased  heating  by 
particle  bombardment  via  increased  charged  particle  number 
densities  at  the  sheath  edge.  Radiative  transfer  in  the  MPD 
plasma  is  a  topic  for  detailed  study.  Although  these  two  de¬ 
ficiencies  can  affect  the  quantitative  predictions  of  the  surface 
temperature,  they  will  not  affect  the  existence  of  the  thermal 
runaway  mechanism  and  unstable  operating  regimes  found 
here.  The  theory  presented  in  this  paper  together  with  models 
of  the  flowing  plasma” 31 32  provide  the  designer  with  a  num¬ 
ber  of  tools.  In  addition  to  the  prediction  of  erosion  rates 
under  a  variety  of  possible  operating  conditions,  the  proper 
cathode  length,  diameter,  material,  external  cooling  condi¬ 
tions.  and  choice  of  propellant  can  be  systematically  evalu¬ 
ated. 

Experimental  verification  of  the  predictions  of  the  theory- 
presented  in  this  paper  requires  simultaneous  measurement 
of  the  plasma  current  density,  plasma  charged  particle  number 
density  near  the  electrode,  amount  of  external  cooling,  and 
the  local  surface  temperature.  Although  temperature  mea¬ 
surements  in  low-power  steady-state  MPD  thrusters  have  been 
made.'1  lack  of  knowledge  with  regard  to  the  local  current 
and  number  densities  makes  comparison  between  this  theory 
and  existing  measurements  unreliable.  The  measured  steady- 
state  cathode  surface  temperatures  in  a  subscale  device  for 
thoriated  tungsten  range  from  2087  to  2281  K  for  a  range  of, 
total  currents  from  500  to  1200  A  and  mass  flows  from  46.4 
to  150.8  mg/s. 34  This  is  comparable  to  the  parameter  range 
considered  in  this  paper.  Despite  the  difficulties  encountered 
in  nonintrusively  measuring  plasma  number  densities  and  cur¬ 
rent  densities,  efforts  to  measure  these  variables  along  with 
the  local  surface  temperature  would  be  desirable  and  inval¬ 
uable.  Such  local  surface  temperature  measurements,  if  based 
on  optical  pvrometry.  must  also  account  for  reflection  from 
the  radiating  MPD  plasma.  As  a  final  point,  it  must  be  men¬ 
tioned  that  scanning  electron  micrographs  of  used  cathodes 
always  display  evidences  of  craters  from  spots  and  microspots 
in  some  regions.  These  spots  can  be  formed  easily  during  the 
startup  transient  or  the  shutdown  transient  of  the  arc.  There¬ 
fore.  future  experiments  aimed  at  understanding  the  diffuse 
mode  limits  must  exercise  great  caution  to  make  sure  that 
spots  are  not  formed  during  arc  initiation  or  shutdown. 
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Nomenclature 

A  =  thermionic  emission  coefficient 
d  =  anode  thickness 

e  =  electronic  charge 

J  =  total  current 

j  =  current  density  passing  through  the  anode 
j,  =  electron  current  density  from  the  plasma  (Its  sense  is 
positive  away  from  the  anode.) 

Je  =  current  density  of  thermionically  emitted  electrons 
(Its  sense  is  positive  into  the  anode.) 
y®  =  average  net  current  density  through  the  plasma,  taken 

positive  away  from  the  anode 
k  -  Boltzmann’s  constant 

L  =  anode  length 

me  =  electron  mass 

ne  =  electron  number  density  at  the  plasma-sheath 
boundary 

0(  )  =  order  of  magnitude  of  the  quantity  in  parenthesis 
<7o  =  heat  flux  on  the  anode  inner  surface  (Its  sense  is 

positive  into  the  anode  surface.) 

T  =  temperature 

T0  =  anode  inner  surface  temperature 

Toe  =  critical  value  of  T0  where  dj„/dT0  =  0 
Td  =  anode  outer  surface  temperature 

Tr  =  plasma  electron  temperature 

VA  =  anode  sheath  voltage  drop,  defined  as  the  potential  at 
the  plasma-sheath  edge  minus  the  anode  potential 
W  =  anode  width 

x  =  outward  coordinate  taken  to  be  0  at  the  anode  inner 

surface,  and  d  at  the  anode  outer  surface 
e  =  emissivity  of  outer  anode  surface 
X  =  thermal  conductivity  of  anode  material 
a sb  =  Stefan-Boltzmanr.  constanr 

a  =  electrical  conductivity  of  anode  material 

4>a  =  anode  material  work  function 


I.  Introduction 

N  important  consideration  in  the  use  of 
dynamic  (MPD)  thrusters  operating  at 


magnetoplasma- 
steady  state  for 
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space  missions  is  electrode  erosion.  Erosion  in  MPD  thrusters 
has  been  studied  experimentally  and  has  focused  mainly  on  the 
cathode.13  There  are  also  ongoing  experiments  in  cathode 
erosion.4  However,  in  this  article,  we  present  a  simple  analysis 
of  the  energy  balance  on  the  anode  in  an  effort  to  predict 
anode  surface  temperatures.  This  can  subsequently  be  used  to 
predict  erosion  rates  by  evaporation.  As  will  be  shown,  even 
such  a  simple  model  reveals  subtle  physical  phenomena. 

Electrode  erosion  is  connected  with  the  modes  of  current 
conduction  through  the  electrode  surface.  Two  modes  are 
known  to  exist:  spot  and  diffuse.  The  diffuse  mode  at  subon¬ 
set  conditions  is  the  focus  of  this  article.  Vainberg  et  al.5  have 
considered  anode  behavior  at  onset  conditions  and  beyond. 
Their  explanation  of  anode  melting  rests  on  the  anode  sheath 
reversal  mechanism,  which  has  also  been  observed  by  Hugel.6 
In  constrast  to  these  earlier  works,  it  is  shown  in  this  article 
that  a  thermal  runaway  may  occur  well  before  the  sheath 
reversal.  This  work  supports  earlier  conjectures  that  local 
surface  melting  of  the  anode  precedes  spot  formation.7  Two 
operating  modes  are  predicted  by  the  theory  presented  here. 
One  of  these  yields  a  stable  steady-state  temperature  for  the 
ariOv’",  whereas  the  other  results  in  a  thermal  runaway  in 
which  the  anode  regeneratively  heats  itself  until  it  melts.  The 
total  current,  anode  geometry,  and  material  work  function  are 
shown  to  strongly  influence  the  steady-state  anode  tempera¬ 
ture  for  the  stable  operating  regime. 

In  the  following  section,  the  governing  equations  describing 
anode  heat  transfer  will  be  derived.  The  solutions  to  these 
equations  under  some  conditions  of  interest  are  given  and 
discussed  in  Sec.  Ill,  followed  by  the  summary  and  conclu¬ 
sions  in  Sec.  IV. 


II.  Anode  Energy  Balance 

Consider  the  hollow  cylindrical  anode  geometry  of  the 
MPD  thruster  modeled  as  a  long  thin  slab  of  length  L ,  width 
W,  and  thickness  d  (shown  in  Fig.  1).  At  steady  state,  the 
energy  balance  gives 


dx2~  a 


(1) 


Equation  (1)  is  subjected  to  the  following  boundary  condi¬ 
tions: 


§ 

II 

o 

M 

* 

ys 

1 

’  (2) 

X  1  X  =  d  ~  «*SB  Ti 

(3) 

For  the  case  of  constant  properties,  the  system  of  Eqs.  (1-3) 
may  be  readily  integrated  to  give 


j2d 

toSBTd  =  —  +  (Jo 


(4) 


and 


\Td 


fd2 

2a 


—Qo/d  +  \7"o 


(5) 
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Fig.  1  The  hollow  cylindrical  anode  is  shown  in  this  idealization  as  a 
long  thin  slab  of  length  L ,  width  W,  and  thickness  d. 


Now,  if  q0  is  known  and  j  is  given,  Eqs.  (4)  and  (5)  represent 
two  equations  for  the  two  unknowns  T0  and  Td. 

The  heat  flux  on  the  anode  inner  surface  qo  must  be  found 
by  considering  particle  bombardment  from  the  plasma 
through  the  sheath.  Consider,  therefore,  a  power  balance  per 
unit  area  on  the  anode  inner  surface.  The  focus  here  is  on 
conditions  prior  to  sheath  reversal;  therefore,  the  anode  is  at 
a  lower  potential  with  respect  to  the  plasma  potential  at  the 
sheath  edge.  The  power  balance,  then,  mainly  consists  of  a 
balance  between  electron  bombardment  and  thermionic  emis¬ 
sion  cooling.  It  can  be  shown  that  under  these  operating 
conditions  evaporative  cooling  is  a  relatively  small  effect.  Ion 
bombardment  may  be  easily  included  in  our  analysis,  but  we 
neglect  it  for  simplicity.  Although  this  is  expected  to  influence 
final  results  quantitatively,  the  conclusions  regarding  the 
mechanisms  responsible  for  a  thermal  runaway  will  not  be 
altered.  We  then  have 

q<1=jA<t>A+2kT'/e)-ATi 


ANODE  INSIDE  TEMPERATURE  (K) 

Fig.  2  A  typical  variation  of  the  heat  defect  (defined  as  the  net  heat 
flux  into  the  anode  minus  the  net  heat  flux  out)  is  shown  here  vs  the 
anode  inside  surface  temperature  7g. 

where  q0  is  given  by  Eq.  (9).  Given  d,  L,  W,  T„  and  7,  the 
axial  current  density  throught  the  anode  ranges  up  to  j-J/ 
Wd,  and  the  average  current  density  through  the  plasma  at  the 
anode  surface  is  jm=J/WL.  Although  this  simplification  is 
unnecessary  (since  j  and  j m  can  be  treated  as  separate  point- 
wise  parameters),  it  allows  the  thermal  steady  state  to  be 
related  directly  to  the  total  current  (a  more  tangible  quantity) 
rather  than  the  current  densities.  With  these  quantities,  Eq. 
(10)  may  be  solved  for  7"0.  Once  T0  is  determined,  Td  may  be 
readily  obtained  from  Eqs.  (4)  or  (5).  Equations  (9)  and  (10) 
have  been  solved  in  this  manner  for  representative  conditions 
in  the  MPD  thruster.  These  results  are  discussed  in  the  follow¬ 
ing  section. 


x  (<*>„  +  2k To/e)  exp  {-e<t>A/kT0\  (6) 

The  electron  current  density  from  the  plasma  je  can  be  deter¬ 
mined  from  overall  current  conservation 

je  =/»  +  A  To  exp(  ~e<t>A/kT0\  (7) 

Inclusion  of  ion  bombardment  would  involve  additional  terms 
on  the  right  hand  side  of  Eqs.  (6)  and  (7).  The  anode  sheath 
drop  may  be  determined  from  Eq.  (7)  to  be 

V*  =  HVm+jeVjr]  (8) 

where  jr=‘ene(kTt/2*me)'A  and  jE  =  A T$  exp\-e<bA/kT0\. 
The  anode  is  at  a  potential  of  —  VA  with  respect  to  the  plasma- 
sheath  edge  (taken  as  the  V  =  0  datum)  in  the  regime  under 
consideration.  Combining  Eqs.  (6)  and  (7)  gives 


.  /  2 kT\  lAkT^T'-To) 

do=J.y>A  +— -J  + - - - expl-e<t>A/kT0l 

(9) 

Finally,  combining  Eqs.  (4)  and  (5),  we  get  an  implicit  equa¬ 
tion  for  T0 

(10) 


III.  Analysis  and  Results 

In  the  previous  section,  it  was  shown  that  an  energy  balance 
on  the  anode  at  steady  state  yields  an  implicit  equation  for  the 
anode  inside  the  surface  temperature  T0.  This  equation  [Eq. 
(10)1  is  nonlinear  in  T0  and  must  be  solved  numerically.  This 
section  will  focus  on  the  solution  of  this  implicit  equation, 
resulting  in  the  discovery  of  a  thermal  runaway  mode  and 
culminating  in  a  discussion  of  various  operating  limits. 

For  illustrative  purposes,  let  us  consider  a  tungsten  anode  of 
width  30  cm  and  length  10  cm.  Let  the  electron  temperature  T, 
be  Fixed  at  15,000  K.  The  typical  variation  of  the  implicit 
function  F(T0)  vs  T0,  given  by  Eq.  (10),  is  displayed  in  Fig.  2. 
The  function  F(T0)  represents  the  net  heat  into  the  anode  per 
unit  area  per  unit  time.  It  is  immediately  evident  from  Fig.  2 
that  two  steady-state  solutions  exist.  Without  performing  a 
stability  analysis,  one  can  conclude  on  the  basis  of  physical 
reasoning  that  the  lower  temperature  is  a  stable  root  and  the 
higher  one  is  an  unstable  root.  Consider  the  smaller  tempera¬ 
ture  First.  Any  perturbation  to  the  left  of  this  root  (i.e.,  on  the 
lower  temperature  side)  results  in  a  net  heat  flux  into  the 
anode.  To  counter  this,  the  surface  temperature  must  increase 
in  order  to  maintain  a  steady  state.  Similarly,  any  small  per¬ 
turbation  to  the  right  (i.e.,  on  the  higher  temperature  side) 
results  in  a  net  heat  flux  out  of  the  anode.  To  counter  this,  the 
surface  temperature  must  decrease.  Thus,  we  see  that  the 
tendency  of  the  thermal  response  is  to  drive  the  temperature 
toward  this  root.  Using  this  physical  argument  on  the  second 
root  (i.e.,  higher  temperature),  we  may  conclude  that  it  is 
unstable.  If  the  surface  temperature  increases  just  above  the 
value  given  by  the  second  root,  the  steady-state  energy  balance 
indicates  a  further  increase  in  temperature  in  order  to  main¬ 
tain  a  steady  state.  This  is  because  the  only  possible  way  of 
losing  energy  is  through  radiation  and  by  thermionic  emission. 
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both  of  which  increase  with  increasing  temperature.  This  re¬ 
sults  in  a  thermal  runaway  with  eventual  melting  of  the  anode. 

A  physical  interpretation  of  the  two  aforementioned  solu¬ 
tions  will  now  be  discussed.  Examination  of  Eq.  (8)  indicates 
that  the  anode  sheath  drop  decreases  as  the  net  plasma  current 
density  increases.  This  results  in  increased  electron  bombard¬ 
ment,  as  can  be  seen  from  Eqs.  (6)  and  (7),  causing  heating  of 
the  anode  surface.  Two  possible  scenarios  may  ensue.  The 
resulting  increase  in  the  surface  temperature  leads  to  an  in¬ 
crease  in  the  emission  current  density,  which  further  lowers 
the  anode  sheath  voltage  drop.  This  positive  feedback  may 
repeat  itself  until  thermal  runaway  leads  to  local  melting  of  the 
surface.  The  second  scenario  is  stable,  steady  operation  arising 
from  the  fact  that  sufficient  cooling  of  the  anode  prevents  the 
thermal  runaway  by  excessive  electron  bombardment.  The  two 
solutions  found  to  Eq.  (10)  represent  these  two  situations. 

The  anode  surface  temperature  70  for  stable  operation  is 
found  to  be  strongly  dependent  on  the  total  current,  anode 
geometry,  and  the  anode  work  function.  Figure  3  shows  the 
anode  thickness  for  various  total  currents,  and  Fig.  4  shows 
the  anode  inside  surface  temperature  vs  total  current  for  vari¬ 
ous  anode  thicknesses.  From  these  results,  it  is  evident  that  for 
a  given  total  current,  there  exists  an  optimum  thickness  for 
achieving  a  minimum  anode  temperature.  An  interesting  fea¬ 
ture  of  the  pair  of  solutions  found  to  the  steady-state  heat 
balance  equation  is  that  the  stable  and  unstable  roots  start 
moving  toward  each  other  as  the  current  is  increased.  This  con¬ 
tinues  until  a  critical  value  of  the  current  is  reached,  beyond 
which  no  steady  solution  can  be  found  below  the  melting  point 
of  the  material.  Also,  for  a  fixed  total  current,  these  roots 
approach  each  other  as  the  thickness  is  increased.  It  is  possible 
that  the  analytical  behavior  of  this  unstable  thermal  runaway 
point  will  yield  insight  into  electrode  material  breakdown. 

Some  interesting  limits  may  be  readily  obtained  by  analyz¬ 
ing  the  equations  presented  in  Sec.  II.  For  instance,  from  Eq. 
(5)  it  is  clear  that 

j2d2/2o<kT0  or  j  <jc  =(2koTa)'/'/d  (11) 


prevent  thermal  runaway.  This  thermal  runaway  can  occur 
under  the  steady,  diffuse  mode  operation  of  the  anode. 

An  important  stability  criterion  may  be  derived  by  requiring 
that  6F/6T0<0  for  stable  diffuse  operation.  Using  Eq.  7,  we 
find 


For  VA  below  the  value  given  by  the  right-hand  side  of  Eq. 
(12),  excessive  electron  bombardment  will  cause  a  thermal 
runaway.  Thus,  local  surface  melting  can  occur  prior  to  sheath 
reversal. 

A  limit  may  also  be  found  for  the  plasma  current  density jm. 
Considering  thin  anodes  (i.e.,  7*  T0  throughout  the  anode), 
neglecting  ohmic  heating  and  any  external  cooling  except  for 
radiation,  we  may  write  the  energy  balance  at  steady  state  as 

ease**  +^y^(7V  -  7b)  exp(  —  e<t>A/kT0) 


Using  the  fact  that  Te>T0,  differentiating  with  respect  to  7b, 
and  setting  djm/dT0=Q,  we  obtain 


(J oe)max 


[1  —(4kToc/e<t>A)] 
[<t>A  +  (2k  7,/e)] 


ettSB7br4 


(14) 


where  is  the  extremum  plasma  current  density,  and  To, 

is  the  critical  value  of  70  where  dj^/dTo-0.  Additional  heat 
transfer  to  an  external  coolant  results  only  in  a  slight  modifi¬ 
cation  of  Eq.  (14): 


Ucd)  max 


[l-(4k70C/e^)l  r  4  [l-jkTpc/eM] 
[<t>A+(2kTe/e))  f°SB  *  l<t>A  +  (2k7,/e)] 


hToc 


which  is  very  similar  (within  a  proportionality  constant)  to  the 
thermal  runaway  condition  obtained  by  Hantzsche8  for 
cathode  spots.  However,  it  must  be  pointed  out  that  Eq.  (11) 
is  obtained  here  for  an  anode  of  constant  electrical  conductiv¬ 
ity,  operating  in  the  diffuse  mode.  This  sharply  differs  from 
Hantzsche’s  consideration  of  a  cathode  spot  with  an  electrical 
conductivity  that  varies  according  to  the  Wiedemann-Franz 
law  (i.e.,  aocl/7).  Furthermore,  in  constrast  to  Hantzsche’s 
work,  surface  cooling  due  to  thermionic  emission  does  not 


(15) 

where  h  is  the  heat  transfer  coefficient  between  the  anode  (at 
To)  and  an  external  coolant  (at  7C  ■<  To)-  For  7,  =  20,000  K, 
Eq.  (14)  approximately  yields  47.6  A/cm2  for  a  tungsten  an¬ 
ode  (<t>A  =4.52  V)  and  5.4  A/cm2  for  a  thoriated  tungsten 
anode  ( <t>A  =2.63  V).  The  lower  work  function  yields  a  lower 
maximum  plasma  current  density  because  of  a  lower  value  of 
Toe-  Since  these  limiting  current  densities  are  far  lower  than 


O  12  3  -V  S  3 

ANODE  THirKNESS  (cm) 

Fig.  3  The  anode  inside  surface  temperature  is  plotted  here  vs  anode 
thickness,  for  various  total  currents. 


Fig.  4  The  anode  inside  surface  temperature  is  shown  here  vs  total 
current,  for  various  anode  thicknesses. 
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those  required  for  a  steady  MPD  discharge,  it  is  clear  that  the 
anode  must  be  cooled  externally.  This  is  done  in  practice.1,4 

IV.  Summary  and  Conclusions 

An  anode  energy  balance  has  been  performed  that  includes 
the  effects  of  the  anode  sheath.  Results  show  that  under  many 
conditions,  two  steady-state  solutions  can  be  found.  One  of 
these  corresponds  to  a  stable  operating  point,  the  other  to  a 
thermal  runaway.  The  stable  root  that  gives  the  anode  inside 
surface  temperature  at  steady  state  was  found  to  be  strongly 
dependent  on  the  total  current,  anode  geometry,  and  the  mate¬ 
rial  work  function.  A  stability  condition  in  the  form  of  a 
minimum  anode  sheath  voltage  drop  [Eq.  (12)]  has  also  been 
given. 

Several  conclusions  can  be  derived  from  the  anode  thermal 
analysis: 

1)  There  exists,  for  a  given  current  density,  an  optimum 
thickness  for  achieving  a  minimum  anode  temperature. 

2)  Under  many  conditions,  there  exists  a  pair  of  steady-state 
solutions  to  the  heat  balance,  one  of  which  is  a  stable  operat¬ 
ing  point,  the  other  of  which  is  an  unstable  thermal  runaway 
point. 

3)  There  exists  a  maximum  steady  operating  current  density 
for  a  given  electrode  geometry  based  on  melting  temperature 
and  cooling  considerations. 

The  thermal  runaway  mode  discovered  for  the  anode  may 
explain  anode  spots  under  MPD  conditions.  Although  some 
limits  have  been  found,  anode  heat  transfer  effects  need  to  be 
explored  further  in  order  to  learn  more  about  what  controls 
the  limits  of  stable  operation.  Since  external  cooling  is  seen  to 
be  extremely  important,  detailed  experimental  results  quan¬ 


tifying  cooling  rates  for  various  operating  conditions  and  ge¬ 
ometries  would  be  valuable.  Radiative  heat  transfer  from  the 
plasma,  which  has  been  neglected  here,  is  another  important 
area  for  further  research. 
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